2023-02-02 Implementing a retinotopic transform using grid_sample from pyTorch

Laurent Perrinet

2023-02-02 08:59

Implementing a retinotopic transform using `grid_sample` from `pyTorch`¶

The grid_sample transform is a powerful function which allows to transform any input image into a new topology. It is notably used in Spatial Transformer Networks for instance to learn CNN to be invariant to affine transforms. We used it recently in a publication What You See Is What You Transform: Foveated Spatial Transformers as a Bio-Inspired Attention Mechanism by Ghassan Dabane et al.

The use of grid_sample can b etedious and here, we show how to use it to create a log-polar transform of the image and create the following figure:

Retinotopy

A picture (extract from the painting "The Ambassadors" by Hans Holbein the Younger can be represented on a regular grid represented by vertical (red) and horizontal (blue) lines. Retinotopy transforms this grid, and in particular the area representing the fovea (shaded gray) is over-represented. Applied to the original image of the portrait, the image is strongly distorted and represents more finally the parts under the axis of sight (here the mouth).

Elections présidentielles 2022: estimation du transfert de voix

Laurent Perrinet

2022-06-08 21:27

tl;dr : On essaie ici de deviner le transfert des voix entre les choix effectués entre deux scrutins d'un vote (ici les élections présidentielles 2022 en France) par une méthode d'apprentissage automatique.

Afin d'analyser les résultats des élections, par exemple les dernières élections présidentielles de 2022 en France, et de mieux comprendre la dynamique des choix de vote entre les différents groupes de population, il peut être utile d'utiliser des outils d'apprentissage automatique pour inférer des structures à première vue cachées dans la masse des données. En particulier, inspiré par cet article du Monde, on peut se poser la question de savoir si on peut extraire depuis les données brutes des élections une estimation du report de voix entre les choix de vote au premier tour et ceux qui sont effectués au deuxième tour.

Pour cela, parmi les outils mathématiques de l'apprentissage automatique, nous allons utiliser des probabilités. Cette théorie va nous permettre d'exprimer le fait que les résultats tels qu'ils sont obtenus peuvent présenter une variabilité mais que celle-ci réelle résulte de préférences de chaque individu dans la population votante. En particulier, on peut considérer que chaque individu va avoir une préférence, graduée entre $0=0\%$ et $1=100\%$ pour chacun des choix (candidats, nul, blanc, abstention) au premier et second tour. Ainsi, les votes effectués vont correspondre à la réalisation de ces préférences.

Bien sûr, le vote reste secret et on n'a pas accès au vote de chaque individu et encore moins à ses préférences. Mais comme chaque bureau de vote présente des variabilités liées au contexte local et qui fait que cette population locale a une préférence pour certains choix plutôt que d'autres, on peut considérer chaque bureau de vote comme une population individuelle pour lequel nous allons essayer de prédire les résultats du vote au deuxième tour. En exploitant les irrégularités locales, bureau de vote par bureau de vote, nous allons pouvoir prédire (le mieux possible) le report des votes individuel (de chaque individu tel qu'il passe d'un vote à un autre, par exemple de "Hidalgo" à "Macron"). Nous allons en particulier montrer une prédiction très correcte des données du second tour à partir de ceux du premier, montrant la plausibilité d'une telle hypothèse :

Prédiction du transfert des voix

C'est à ma connaissance une contribution originale (jusqu'à ce qu'une bonne âme veuille bien me donner un lien vers une méthode existante similaire qui me permette de correctement la citer...) que nous allons exploiter ici. Cette prédiction, si elle est efficace (et on va montrer qu'elle est en moyenne correctement prédite avec moins de 2 points de pourcentages d'erreur près), peut donner une idée du transfert de vote entre les deux tours qui a lieu en fonction des préférences des votes de chaque individu.

Nous allons dans la suite montrer comment on peut estimer le pourcentage de chances d'exprimer une voix pour un candidat ou pour l'autre en fonction du choix qu'on a exprimé au premier tour:

Transfert des voix

Comme on le verra plus bas, ce tableau montre des tendances claires, par exemple que si on a voté "Macron", "Jadot", "Hidalgo" ou "Pécresse" au premier tour, alors on va certainement voter "Macron" au deuxième tour. Ces électeurs se montrent particulièrement consensuel et suivent le « pacte républicain » mise en place pour faire un "barrage" au Front National (en suivant le terme consacré). Il montre aussi que si on a voté "Le Pen" ou "Dupont-Aignan" au premier tour alors on va voter Le Pen au deuxième, un clair vote de suivi.

Connaissant les couleurs politiques d'autres candidats du premier tour, on peut être surpris que les électeurs de "Arthaud", "Roussel", "Lassalle" ou "Poutou" ont majoritairement choisi "Le Pen" au deuxième tour, signifiant alors un rejet du candidat Macron. Les électeurs de Zemmour sont aussi partagés, signifiant un rejet des deux alternatives. Ce résultat est à prendre avec des pincettes car ces derniers candidats ont obtenu moins de votes et donc que le processus d'inférence est forcément moins précis car il y a moins de données disponibles.

En résumé, cette analyse donne des tendances en fonction des choix exprimés au premier tour: Transfert des voix qui montre une nette séparation des groupes de vote.

COSYNE reviewer feedback

Laurent Perrinet

2022-02-11 12:35

tl;dr : Crowd-sourcing raw scores for your COSYNE reviewer feedback.

Following that message:

Dear community,

COSYNE is a great conference which plays a pivotal role in our field. If you have submitted an abstract (or several) you have recently received your scores. I am not affiliated to COSYNE - yet willing to contribute in some way: I would like to ask one minute of your time to report the raw scores from your reviewers. I will summarize in a few lines the results in one week time (11/02). The more numerous your feedbacks the better their precision!

Thanks!

As of 2022-02-20, I had received $N = 98$ answers from the google form (out of them, $95$ are valid) out of the $881$ submitted abstracts. In short, the result is that the total score $S$ is simply the linear sum of the scores $s_i$ given by each reviewer $i$ relatively weighted by the confidence levels $\pi_i$ (as stated in the email we received from the chairs):

$$ S = \frac{ \sum_i \pi_i\cdot s_i}{\sum_i \pi_i} $$

Or if you prefer $$ S = \sum_i \frac{\pi_i}{\sum_j \pi_j} \cdot s_i $$

We deduce from that formula that the threshold is close to $6.34$ this year:

2022-02-11_COSYNE-razor

More details in the notebook (or directly in this post) which can also be forked here and interactively modified on binder.

EDIT: On 2022-02-20, I have updated the notebook to account for new answers, I have now received $N = 98$ answers (out of them, $95$ are valid), yet nothing changed qualitatively. On 2022-02-11, I had received $N = 82$ answers from the google form (out of them, $79$ are valid) and the estimated threshold wass close to $6.05$.

Dreamachine

Laurent Perrinet

2022-01-30 11:33

It's at the MIAM (Miam Musée International des Arts Modestes) in Sète, France, that I could for the first time experience really the Dreamachine. It's an optical system which consists of a central light which is periodically occluded by a rotating (cardboard?) cylinder.

The magic of it is that the frequency of occlusion is around $12$ Hz, an important resonant state for sensory system. For the first time, I could really try it out at the MIAM - the important point being to close your lids and rest quiet while looking at the stroboscopic light source. Surprisingly, you see the emergence of "psychedelic patterns" (of course, less than in hippie's movies) yet of the order of the color pattern that may arise in Benham's Disk.

It's difficult to reproduce this pattern on a screen, yet it is still possible to give an impression of it. The goal is here :

to generate a complex visual stimulation flickering on average at $12$ Hz
to project it on a retinotopic space to maximise the "psychelic" effect

retino_alpha

Logistic regression

Laurent Perrinet

2021-12-12 10:32

In a previous notebook, I have shown how to fit a psychometric curve using pyTorch. Here, I would like to review some nice properies of the logistic regression model (WORK IN PROGRESS).

Experimenting with transfer learning for visual categorization

Laurent Perrinet

2021-12-03 16:23

Hi! I am Jean-Nicolas Jérémie and the goal of this notebook is to provide a framework to implement (and experiment with) transfer learning on deep convolutional neuronal network (DCNN). In a nutshell, transfer learning allows to re-use the knowlegde learned on a problem, such as categorizing images from a large dataset, and apply it to a different (yet related) problem, performing the categorization on a smaller dataset. It is a powerful method as it allows to implement complex task de novo quite rapidly (in a few hours) without having to retrain the millions of parameters of a DCNN (which takes days of computations). The basic hypothesis is that it suffices to re-train the last classification layers (the head) while keeping the first layers fixed. Here, these networks teach us also some interesting insights into how living systems may perform such categorization tasks.

Based on our previous work, we will start from a VGG16 network loaded from the torchvision.models library and pre-trained on the Imagenet dataset wich allows to perform label detection on naturals images for $K = 1000$ labels. Our goal here will be to re-train the last fully-Connected layer of the network to perfom the same task but in a sub-set of $K = 10$ labels from the Imagenet dataset.

Moreover, we are going to evaluate different strategies of transfer learning:

VGG General : Substitute the last layer of the pyTorch VGG16 network ($K = 1000$ labels) with a new layer build from a specific subset ($K = 10$ labels).
VGG Linear : Add a new layer build from a specific subset ($K = 10$ labels) after the last Fully-Connected layer of the the pyTorch VGG16 network.
VGG Gray : Same architecture as the VGG General network but trained with grayscale images.
VGG Scale : Same architecture as the VGG General network but trained with images of different size.
VGG Full : Same architecture as the VGG General network but all the layers are trained (otherwise I trained the last Fully-Connected layer).

In this notebook, I will use the pyTorch library for running the networks and the pandas library to collect and display the results. This notebook was done during a master 2 internship at the Neurosciences Institute of Timone (INT) under the supervision of Laurent Perrinet. It is curated in the following github repo.

Find a 8 in a forest of 9

Laurent Perrinet

2021-12-01 21:30

The goal here is to find one non-repeting pattern in an image. First, find the pattern, then second, compute correlation.

Triangulating stars on the night sky

Laurent Perrinet

2021-12-01 10:32

In a previous notebook, I have shown properties of the distribution of stars in the sky. Here, I would like to use the existing database of stars' positions and display them as a triangulation

Modelling Ocean surface waves

Laurent Perrinet

2021-09-11 11:56

Trying to model Ocean surface waves, I have played using a linear superposition of sinusoids and also on the fact that phase speed (following a wave's crest) is twice as fast as group speed (following a group of waves). Finally, I more recently used such generation of a model of the random superposition of waves to model the wigggling lines formed by the refraction (bendings of the trajectory of a ray at the interface between air and water) of light, called caustics...

Observing real-life ocean waves instructed me that if a single wave is well approximated by a sinusoïd, each wave is qualitatively a bit different. Typical for these surface gravity waves are the sharp crests and flat troughs. As a matter of fact, modelling ocean waves is on one side very useful (Ocean dynamics and its impact on climate, modelling tides, tsunamis, diffraction in a bay to predict coastline evolution, ...) but quite demanding despite a well known mathematical model. Starting with the Navier-Stokes equations to an incompressible fluid (water) in a gravitational field leads to Luke's variational principle in certain simplifying conditions. Further simplifications lead to the approximate solution given by Stokes which gives the following shape as the sum of different harmonics:

This seems well enough at the moment and I will capture this shape in this notebook and notably if that applies to a random mixture of such waves...

(WORK IN PROGRESS)

The current situation is that these solutions seem to not fit what is displayed on the wikipedia page and that I do not spot the bug I may have introduced... More work on this is needed and any feedback is welcome...

Density of stars on the surface of the sky

Laurent Perrinet

2021-03-27 14:06

Looking at the night sky, the pattern of stars on the surface of the sky follows a familiar pattern. The Big Dipper, Cassiopeia, the Pleiades, or Orion are popular landmarks in the sky which we can immediatly recognize.

Different civilisations labelled these patterns using names such as constellations in the western world. However, this pattern is often the result of pure chance. Stars from one constellation belong often to remote areas in the universe and they bear this familiarity only because we always saw them as such; the rate at which stars move is much shorter than the lifespan of humanity.

I am curious here to study the density of stars as they appear on the surface of the sky. It is both a simple question yet a complex to formulate. Is there any generic principle that could be used to characterize their distribution? This is my attempt to answer the question

https://astronomy.stackexchange.com/questions/43147/density-of-stars-on-the-surface-of-the-sky

(but also to make the formulation of the question clearer...). This is my answer:

Time lapsing an orchid's flower

Laurent Perrinet

2021-02-21 20:57

The goal here is to realize a time-lapse using a raspberry π and some python code and finally get this:

(This was done over 10 days, with almost each day an irregular session the morning and one in the evening)

Kuramoto model

Laurent Perrinet

2021-01-08 14:42

The goal here is to re-implement the Kuramoto model, following a lecture from Joana Cabral and the code that is provided, but using python instead of matlab.

Generating second-order figures from texture

Laurent Perrinet

2021-01-08 11:54

The goal here is to use the MotionClouds library to generate figures with second-order contours, similar to those used in the P. Roelfsema's group.

Extracting music from the screenshots of a Spotify playlist

Laurent Perrinet

2020-11-24 15:19

In this post I'll try to show how from a screenshot obtained from a software like Spotify you can programmatically extract the tracks of the songs as well as the artists, to finally download them from the Internet.

Fitting COVID data

Laurent Perrinet

2020-10-10 14:13

I propose here a simple method to fit experimental data common to epidemiological spreads, such as the present COVID-19 pandemic, using the inverse gaussian distribution. This follows the general incompregension of my answer to the question Is the COVID-19 pandemic curve a Gaussian curve? on StackOverflow. My initial point is to say that a Gaussian is not adapted as it handles a distribution on real numbers, while such a curve (the variable being number of days) handles numbers on the half line. Inspired by the excellent A Theory of Reaction Time Distributions by Dr Fermin Moscoso del Prado Martin a constructive approach is to propose another distribution, such as the inverse Gaussian distribution.

This notebook develops this same idea on real data and proves numerically how bad the Gaussian fit is compared to the latter. Thinking about the importance of doing a proper inference in such a case, I conclude

Benchmarking CNNs

Laurent Perrinet

2020-09-28 17:25

Hi! I am Jean-Nicolas Jérémie and the goal of this benchmark is to offer a comparison between differents pre-trained image recognition's networks based on the Imagenet dataset wich allows to work on naturals images for $1000$ labels. These different networks tested here are taken from the torchvision.models library : AlexNet, VGG16, MobileNetV2 and ResNet101.

Our use case is to measure the performance of a system which receives a sequence of images and has to make a decision as soon as possible, hence with batch_size=1. Specifically, we wish also to compare different computing architectures such as CPUs, desktop GPUs or other more exotic platform such as the Jetson TX2 (experiment 1). Additionally, we will implement some image transformations as up/down-sampling (experiment 2) or transforming to grayscale (experiment 3) to quantify their influence on the accuracy and computation time of each network.

In this notebook, I will use the Pytorch library for running the networks and the pandas library to collect and display the results. This notebook was done during a master 1 internship at the Neurosciences Institute of Timone (INT) under the supervision of Laurent PERRINET. It is curated in the following github repo.

nesting jupyter runs

Laurent Perrinet

2020-08-09 10:52

Jupyter notebooks are a great way of sharing knowledge in science, art, programming. For instance, in a recent musing, I tried to programmatically determine the color of the sky. This renders as a web page, but is also a piece of runnable code.

As such, they are also great ways to store the knowledge that was acquired at a given time and that could be reusable. This may be considered as bad programming and may have downsides as described in that slides :

Recently, thanks to an answer to a stack overflow question, I found a way to overcome this by detecting if the caall to a notebook is made from the notebook itself or from a parent.

Colors of the sky

Laurent Perrinet

2020-07-04 12:05

Our sensorial environment contains multiple regularities which our brain uses to optimize its representation of the world: objects fall most of the time downwards, the nose is usually in the middle below the eyes, the sky is blue... Concerning this last point, I wish here to illustrate the physical origins of this phenomenon and in particular the range of colors that you may observe in the sky.

Caustic (optics)

Laurent Perrinet

2020-06-19 10:01

Caustics (wikipedia are luminous patterns which are resulting from the superposition of smoothly deviated light rays. It is for instance the heart-shaped pattern in your cup of coffee which is formed as the rays of the sun are reflected on the cup's surface. It is also the wiggly pattern of light curves that you will see on the floor of a pool as the sun's light is refracted at the surface of the water. Here, we simulate that particular physical phenomenon. Simply because they are mesmerizingly beautiful, but also as it is of interest in visual neuroscience. Indeed, it speaks to how images are formed (more on this later), hence how the brain may understand images.

In this post, I will develop a simple formalism to generate such patterns, with the paradoxical result that it is very simple to code yet generates patterns with great complexity, such as:

This is joint work with artist Etienne Rey, in which I especially follow the ideas put forward in the series Turbulence.

Creating an hexagonal grid

Laurent Perrinet

2020-04-16 10:38

A quick note as to how to create an hexagonal grid.

Fitting a psychometric curve using pyTorch

Laurent Perrinet

2020-04-08 16:29

The goal here is to compare methods which fit data with psychometric curves using logistic regression. Indeed, after (long) experiments where for instance you collected sequences of keypresses, it is important to infer at best the parameters of the underlying processes: was the observer biased, or more precise?

While I was forevever using sklearn or lmfit (that is, scipy's minimize) and praised these beautifully crafted methods, I sometimes lacked some flexibility in the definition of the model. This notebook was done in collaboration with Jenna Fradin, master student in the lab.

tl; dr = Do not trust the coefficients extracted by a fit without validating for methodological biases.

One part of flexibility I missed is taking care of the lapse rate, that is the frequency with which you just miss the key. In a psychology experiment, you often see a fast sequence of trials for which you have to make a perceptual deccision, for instance press the Left or Right arrows. Sometimes you know the answer you should have done, but press the wrong eror. This error of distraction is always low (in the order of 5% to 10%) but could potentially change the results of the experiments. This is one of the aspects we will evaluate here.

In this notebook, I define a fitting method using pytorch which fits in a few lines of code :

In [1]:

import torch
from torch.utils.data import TensorDataset, DataLoader

torch.set_default_tensor_type("torch.DoubleTensor")
criterion = torch.nn.BCELoss(reduction="sum")

class LogisticRegressionModel(torch.nn.Module):
    def __init__(self, bias=True, logit0=-2, theta0=0, log_wt=torch.log(0.1*torch.ones(1))):
        super(LogisticRegressionModel, self).__init__()
        #self.linear = torch.nn.Linear(1, 1, bias=bias)
        self.theta0 = torch.nn.Parameter(theta0 * torch.ones(1))
        self.logit0 = torch.nn.Parameter(logit0 * torch.ones(1))
        self.log_wt = torch.nn.Parameter(log_wt * torch.ones(1))

    def forward(self, theta):
        p0 = torch.sigmoid(self.logit0)
        #out = p0 / 2 + (1 - p0) * torch.sigmoid(self.linear(theta))
        out = p0 / 2 + (1 - p0) * torch.sigmoid((theta-self.theta0 )/torch.exp(self.log_wt))
        return out

learning_rate = 0.005
beta1, beta2 = 0.9, 0.999
betas = (beta1, beta2)
num_epochs = 2 ** 9 + 1
batch_size = 256
amsgrad = False # gives similar results
amsgrad = True  # gives similar results

def fit_data(
    theta,
    y,
    learning_rate=learning_rate,
    batch_size=batch_size,  # gamma=gamma,
    num_epochs=num_epochs,
    betas=betas,
    verbose=False, **kwargs
):

    Theta, labels = torch.Tensor(theta[:, None]), torch.Tensor(y[:, None])
    loader = DataLoader(
        TensorDataset(Theta, labels), batch_size=batch_size, shuffle=True
    )

    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

    logistic_model = LogisticRegressionModel()
    logistic_model = logistic_model.to(device)
    logistic_model.train()
    optimizer = torch.optim.Adam(
        logistic_model.parameters(), lr=learning_rate, betas=betas, amsgrad=amsgrad
    )
    for epoch in range(int(num_epochs)):
        logistic_model.train()
        losses = []
        for Theta_, labels_ in loader:
            Theta_, labels_ = Theta_.to(device), labels_.to(device)
            outputs = logistic_model(Theta_)
            loss = criterion(outputs, labels_)

            optimizer.zero_grad()
            loss.backward()
            optimizer.step()
            losses.append(loss.item())

        if verbose and (epoch % (num_epochs // 32) == 0):
            print(f"Iteration: {epoch} - Loss: {np.sum(losses)/len(theta):.5f}")

    logistic_model.eval()
    Theta, labels = torch.Tensor(theta[:, None]), torch.Tensor(y[:, None])
    outputs = logistic_model(Theta)
    loss = criterion(outputs, labels).item() / len(theta)
    return logistic_model, loss

and run a series of tests to compare both methods.

Changing the global phase of a Motion Cloud

Laurent Perrinet

2020-01-08 21:43

Motion Clouds were defined in the origin to define parameterized moving textures. In that other post, we defined a simple code to generate static images using a simple code. Can we generate a series of images while changing the phase globally?

2020-01-08-MC_phase

Role of gamma correction in Sparse coding

Laurent Perrinet

2019-11-28 12:46

I have previously shown a python implementation which allows for the extraction a sparse set of edges from an image. We were using the raw luminance as the input to the algorithm. What happens if you use gamma correction?

albert on gamma

Results : for this particular image, we checked that using the luminance ($\gamma \approx 1$) is the correct choice. The outcome is that gamma correction may improve coding, but only slightly. In the figure below, we plot as a function of gamma the final energy of the error and the perceptually relevant measure of structural simailarity :

albert on gamma

This may not be the case for other types of images which would justify an image-by-image local gain control.
more information on sparse coding is exposed in the following book chapter (see also https://laurentperrinet.github.io/publication/perrinet-15-bicv ):

@inbook{Perrinet15bicv,
    title = {Sparse models},
    author = {Perrinet, Laurent U.},
    booktitle = {Biologically-inspired Computer Vision},
    chapter = {13},
    editor = {Keil, Matthias and Crist\'{o}bal, Gabriel and Perrinet, Laurent U.},
    publisher = {Wiley, New-York},
    year = {2015}
}

Converting a bunch of files in a few lines

Laurent Perrinet

2019-11-19 14:40

You may have a bunch of files that you want to convert from one format to another: images, videos, music, text, ... How do you convert them while using ZSH as your shell language in a single line?

I will take the example of music files which I wish totransform from FLAC to OPUS.

Neurostories: creating videos of the flash-lag effect

Laurent Perrinet

2019-09-26 16:29

This page recollects some variations on the flash-lag effect... and the way to easily and programmmatically generate those movies of the illusion.

Video abstract for "An Adaptive Homeostatic Algorithm for the Unsupervised Learning of Visual Features"

Laurent Perrinet

2019-09-11 10:06

This video shows the results of unsupervised learning with different type of kezrnel normalization. This is to illustrate the results obtained in this paper on the An Adaptive Homeostatic Algorithm for the Unsupervised Learning of Visual Features which is now in press.

Mainen & Sejnowski, 1995

Laurent Perrinet

2019-05-28 10:06

Le but ici de cette première tache est de créer un "raster plot" qui montre la reproducibilité d'un train de spike avec des répétitions du même stimulus. En particulier, nous allons essayer de répliquer la figure 1 de Mainen & Sejnowski (1995).

Ce notebook a été élaboré lors d'un TP dans le cadre du Master 1 de sciences cognitives de l'Université d'Aix-Marseille.

Generating textures with different complexities

Laurent Perrinet

2019-01-30 10:06

This video shows different MotionClouds with different complexities, from a crystal-like Grating to textures with an incresing span of spatial frequencies (resp "Mc Narrow" and "MC Broad"). This is to illustrate the different stimuli used in this paper on the chracterization of speed-selectivity in the retina available @ https://www.nature.com/articles/s41598-018-36861-8 .

Origins of the Von Mises distribution

Laurent Perrinet

2018-12-23 20:31

The goal here is to check if the Von Mises distribution is the a priori choice to make when handling polar coordinates.

Sparse coding of large images

Laurent Perrinet

2018-12-13 10:53

In this post we would try to show how one could infer the sparse representation of an image knowing an appropriate generative model for its synthesis. We will start by a linear inversion (pseudo-inverse deconvolution), and then move to a gradient descent algorithm. Finally, we will implement a convolutional version of the iterative shrinkage thresholding algorithm (ISTA) and its fast version, the FISTA.

For computational efficiency, all convolutions will be implemented by a Fast Fourier Tansform, so that a standard convolution will be mathematically exactly similar. We will benchmark this on a realistic image size of $512 \times 512$ giving some timing results on a standard laptop.

Feature vs global motion

Laurent Perrinet

2018-11-29 10:53

As we see a visual scene, there is contribution of the motion of each of the objects that constitute the visual scene into detecting its global motion. In particular, it is debatable to know which weight individual features, such as small objects in the foreground, have into this computation compared to a dense texture-like stimulus, as that of the background for instance.

Here, we design a a stimulus where we control independently these two aspects of motions to titrate their relative contribution to the detection of motion.

Can you spot the motion ? Is it more going to the upper left or to the upper right?

(For a more controlled test, imagine you fixate on the center of the movie.)

Embedding a trajectory in noise

Laurent Perrinet

2018-11-13 11:43

Motion detection has many functions, one of which is the potential localization of the biomimetic camouflage seen in this video. Can we test this in the lab by using synthetic textures such as MotionClouds?

However, by construction, MotionClouds have no spatial structure and it seems interesting to consider more complex trajectories. Following a previous post, we design a trajectory embedded in noise.

Can you spot the motion ? from left to right or reversed ?

(For a more controlled test, imagine you fixate on the top left corner of each movie.)

(Upper row is coherent, lower row uncoherent / Left is -> Right column is <-)

Rainbow effect

Laurent Perrinet

2018-11-11 13:08

When I see a rainbow, I perceive the luminance inside the arc to be brighter than outside the arc. Is this effect percpetual (inside our head) or physical (inside each droplet in the sky). So, this is a simple notebook to show off how to synthesize the image of a rainbow on a realistic sky. TL;DR: there must be a physical reason for it.

Outline: The rainbow is a set of colors over a gradient of hues, masked for certain ones. The sky will be a gradient over blueish colors.

Statistics of the natural input to a ring model

Laurent Perrinet

2018-11-05 15:41

What does the input to a population of neurons in the primary visual cortex look like? In this post, we will try to have a feeling of the structure and statistics of the natural input to such a "ring" model.

This notebook explores this question using a retina-like temporal filtering and oriented Gabor-like filters. It produces this polar plot of the instantaneous energy in the different orientations for a natural movie :

One observes different striking features in the structure of this input to populations of V1 neurons:

input is sparse: often, a few orientations dominate - the shape of these components (bandwidth) seem to be similar,
there are many "switches": at some moments, the representations flips to another. This is due to cuts in the movie (changes from one scene to the other for instance). In a more realistic setting where we would add eye movements, these switches should also happen during saccades (but is there any knowledge of the occurence of the switch by the visual system?),
between switches, there is some coherence in amplitude (a component will slowly change its energy) but also in time (a component is more likely to have a ghradually changing oriantation, for instance when the scene rotates).

This structure is specific to the structure of natural images and to the way they transform (translations, rotations, zooms due to the motion and deformation of visual objects). This is certainly incorporated as a "prior" information in the structure of the visual cortex. As to know how and where this is implemented is an open scientific question.

This is joint work with Hugo Ladret.

Extending datasets in pyTorch

Laurent Perrinet

2018-09-07 12:27

PyTorch is a great library for machine learning. You can in a few lines of codes retrieve a dataset, define your model, add a cost function and then train your model. It's quite magic to copy and paste code from the internet and get the LeNet network working in a few seconds to achieve more than 98% accuracy.

However, it can be tedious sometimes to extend existing objects and here, I will manipulate some ways to define the right dataset for your application. In particular I will modify the call to a standard dataset (MNIST) to place the characters at random places in a large image.

Generating an unique seed for a given filename

Laurent Perrinet

2018-06-13 13:02

When creating large simulations, you may sometimes create unique identifiers for each of it. This is useful to cache intermediate results for instance. This is the main function of hashes. We will here create a simple one-liner function to generate one.

Testing more complex trajectories

Laurent Perrinet

2018-01-16 10:53

MotionClouds may be considered as a control stimulus - it seems more interesting to consider more complex trajectories.

Predictive coding of variable motion

Laurent Perrinet

2017-12-21 07:58

In some recent modeling work:

Laurent Perrinet, Guillaume S. Masson. Motion-based prediction is sufficient to solve the aperture problem. Neural Computation, 24(10):2726--50, 2012 https://laurentperrinet.github.io/publication/perrinet-12-pred

we study the role of transport in modifying our perception of motion. Here, we test what happens when we change the amount of noise in the stimulus.

In this script the predictive coding is done using the MotionParticles package and for a motion texture within a disk aperture.

accessing the data from a pupil recording

Laurent Perrinet

2017-12-13 13:22

I am experimenting with the pupil eyetracker and could set it up (almost) smoothly on a macOS. There is an excellent documentation, and my first goal was to just record raw data and extract eye position.

In [1]:

from IPython.display import HTML
HTML('<center><video controls autoplay loop src="https://laurentperrinet.github.io/sciblog/files/2017-12-13_pupil%20test_480.mp4" width=61.8%/></center>')

Out[1]:

This video shows the world view (cranio-centric, from a head-mounted camera fixed on the frame) with overlaid the position of the (right) eye while I am configuring a text box. You see the eye fixating on the screen then jumping somewhere else on the screen (saccades) or on the keyboard / hands. Note that the screen itself shows the world view, such that this generates an self-reccurrent pattern.

For this, I could use the capture script and I will demonstrate here how to extract the raw data in a few lines of python code.

Computing sparseness of natural images with retina-like RFs

Laurent Perrinet

2017-11-21 11:56

In [1]:

%load_ext autoreload
%autoreload 2

SparseEdges : computing sparseness of natural images with retina-like RFs¶

Let's compute the "edges" produced with symmetrical filters.

Initialization¶

defining framework

In [2]:

from __future__ import division, print_function
pltimport numpy as np
np.set_printoptions(precision=2, suppress=True)

cluster = False

experiment = 'retina-sparseness'
name_database = 'serre07_distractors'
#parameter_file = '/Users/laurentperrinet/pool/science/BICV/SparseEdges/default_param.py'
parameter_file = 'https://raw.githubusercontent.com/bicv/SparseEdges/master/default_param.py'
#lena_file = '/Users/laurentperrinet/pool/science/BICV/SparseEdges//database/lena256.png'
lena_file = 'https://raw.githubusercontent.com/bicv/SparseEdges/master/database/lena256.png'
lena_file = '../../BICV/SparseEdges/database/lena256.png'
N_image = 100
N = 2**11
B_theta = np.inf
do_linear = False

In [3]:

from SparseEdges import SparseEdges
mp = SparseEdges(parameter_file)
mp.pe.N_X, mp.pe.N_Y = 64, 64
mp.pe.figpath, mp.pe.formats, mp.pe.dpi = 'figures', ['png', 'pdf', 'jpg'], 450
mp.init()
print ('Range of spatial frequencies: ', mp.sf_0)

Range of spatial frequencies:  [ 0.62  0.38  0.24  0.15  0.09  0.06  0.03  0.02]

In [4]:

mp.pe

Out[4]:

{'B_sf': 0.4,
 'B_theta': 0.17453277777777776,
 'C_range_begin': -5,
 'C_range_end': 10.0,
 'MP_alpha': 0.7,
 'MP_do_mask': True,
 'MP_rho': None,
 'N': 2048,
 'N_Dtheta': 24,
 'N_X': 64,
 'N_Y': 64,
 'N_image': None,
 'N_phi': 12,
 'N_r': 6,
 'N_scale': 5,
 'N_svm_cv': 50,
 'N_svm_grid': 32,
 'base_levels': 1.618,
 'd_max': 2.0,
 'd_min': 0.5,
 'd_width': 45.0,
 'datapath': 'database',
 'dip_B_psi': 0.1,
 'dip_B_theta': 1.0,
 'dip_epsilon': 0.5,
 'dip_scale': 1.5,
 'dip_w': 0.2,
 'do_edgedir': False,
 'do_mask': True,
 'do_rank': False,
 'do_whitening': True,
 'dpi': 450,
 'edge_mask': True,
 'edge_scale_chevrons': 180.0,
 'edgefigpath': 'results/edges',
 'edgematpath': 'data_cache/edges',
 'eta_SO': 0.0,
 'figpath': 'figures',
 'figsize': 14.0,
 'figsize_cohist': 8,
 'figsize_hist': 8,
 'formats': ['png', 'pdf', 'jpg'],
 'gamma_range_begin': -14,
 'gamma_range_end': 3,
 'kappa_phase': 0.0,
 'line_width': 1.0,
 'line_width_chevrons': 0.75,
 'loglevel_max': 7,
 'mask_exponent': 3.0,
 'matpath': 'data_cache',
 'multiscale': True,
 'n_theta': 24,
 'noise': 0.33,
 'scale': 0.8,
 'scale_chevrons': 2.5,
 'scale_circle': 0.08,
 'scale_invariant': True,
 'seed': 42,
 'svm_KL_m': 0.34,
 'svm_log': False,
 'svm_max_iter': -1,
 'svm_n_jobs': 1,
 'svm_norm': False,
 'svm_test_size': 0.2,
 'svm_tol': 0.001,
 'verbose': 30,
 'weight_by_distance': True,
 'white_N': 0.07,
 'white_N_0': 0.0,
 'white_alpha': 1.4,
 'white_f_0': 0.4,
 'white_n_learning': 0,
 'white_name_database': 'serre07_distractors',
 'white_recompute': False,
 'white_steepness': 4.0}

In [2]:

import matplotlib
pylab_defaults = { 
    'font.size': 10,
    'xtick.labelsize':'medium',
    'ytick.labelsize':'medium',
    'text.usetex': False,
    'font.family' : 'sans-serif',
    'font.sans-serif' : ['Helvetica'],
    }
matplotlib.rcParams.update(pylab_defaults)

%matplotlib inline
import matplotlib.pyplot as plt
%config InlineBackend.figure_format='retina'
#%config InlineBackend.figure_format = 'svg'
fig_width_pt = 397.48  # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27               # Convert pt to inches
fig_width = fig_width_pt*inches_per_pt  # width in inches
#fig_width = 21
figsize=(fig_width, .618*fig_width)
#figpath, ext = os.path.join(os.getenv('HOME'), 'pool/science/RetinaClouds/2016-05-20_nips'), '.pdf'

Standard edges are oriented, but one may modify that:

In [6]:

sf_0 = .09 # TODO .1 cycle / pixel (Geisler)
params= {'sf_0':sf_0, 'B_sf': mp.pe.B_sf, 'theta':np.pi, 'B_theta': mp.pe.B_theta}
FT_lg = mp.loggabor(mp.pe.N_X/2, mp.pe.N_Y/2, **params)
#(fourier_domain(mp.normalize(np.absolute(FT_lg), center=False))+ image_domain(mp.normalize(mp.invert(FT_lg), center=False)))
fig, a1, a2 = mp.show_FT(FT_lg, axis=True, figsize=(fig_width, fig_width/2))
fig.tight_layout()
mp.savefig(fig, experiment + '_loggabor')

In [7]:

sf_0 = .06 # TODO .1 cycle / pixel (Geisler)
params= {'sf_0':sf_0, 'B_sf': mp.pe.B_sf, 'theta':0., 'B_theta': np.inf}
FT_lg = mp.loggabor(mp.pe.N_X/2, mp.pe.N_Y/2, **params)
fig, a1, a2 = mp.show_FT(FT_lg, axis=True, figsize=(fig_width, fig_width/2))
fig.tight_layout()
mp.savefig(fig, experiment + '_dog')

When defining the framework, one thus needs only one angle:

In [8]:

print ('Range of angles (in degrees): ', mp.theta*180./np.pi)
mp.pe.n_theta = 1
mp.pe.B_theta = np.inf
mp.init()
print ('Range of angles (in degrees): ', mp.theta*180./np.pi)

Range of angles (in degrees):  [-82.5 -75.  -67.5 -60.  -52.5 -45.  -37.5 -30.  -22.5 -15.   -7.5   0.
   7.5  15.   22.5  30.   37.5  45.   52.5  60.   67.5  75.   82.5  90. ]
Range of angles (in degrees):  [ 90.]

In [9]:

print('Final sparseness in the representation = {}'.format(mp.pe.N/mp.oc))
print('Final sparseness in the pyramid = {}'.format(mp.pe.N/(4/3*mp.pe.N_X*mp.pe.N_Y)))

Final sparseness in the representation = 0.0625
Final sparseness in the pyramid = 0.375

one example image¶

In [10]:

mp = SparseEdges(parameter_file)
mp.pe.figpath, mp.pe.formats, mp.pe.dpi = 'figures', ['png', 'pdf', 'jpg'], 450

image = mp.imread(lena_file)
mp.pe.N = N
mp.pe.do_mask = True
mp.pe.n_theta = 1
mp.pe.B_theta = B_theta
mp.pe.line_width = 0
mp.pe.mask_exponent = 4.
mp.init()
image = mp.normalize(image, center=False)
image *= mp.mask
print(image.min(), image.max())
fig, ax = mp.imshow(image, mask=True, norm=False)

-0.950759120364 0.891900615494

In [11]:

name = experiment.replace('sparseness', 'lena')
matname = os.path.join(mp.pe.matpath, name + '.npy')
try:
    edges = np.load(matname)
except:
    edges, C_res = mp.run_mp(image, verbose=False)
    np.save(matname, edges)    

    
matname = os.path.join(mp.pe.matpath, name + '_rec.npy')
try:
    image_rec = np.load(matname)
except:
    image_rec = mp.reconstruct(edges, mask=True)        
    np.save(matname, image_rec)

In [12]:

print(matname)

data_cache/retina-lena_rec.npy

In [13]:

#mp.pe.line_width = 0
fig, a = mp.show_edges(edges, image=mp.dewhitening(image_rec), show_phase=False, mask=True)
mp.savefig(fig, name)

In [14]:

#list_of_number_of_edge =  np.logspace(0, 11, i, base=2)
#list_of_number_of_edge =  4**np.arange(6)
list_of_number_of_edge =  2* 4**np.arange(6)
list_of_number_of_edge =  64* 2**np.arange(6)
print(list_of_number_of_edge)

[  64  128  256  512 1024 2048]

In [15]:

fig, axs = plt.subplots(1, len(list_of_number_of_edge), figsize=(3*fig_width, 3*fig_width/len(list_of_number_of_edge)))
vmax = 1.
image_rec = mp.reconstruct(edges, mask=True)        
vmax = mp.dewhitening(image_rec).max()
for i_ax, number_of_edge in enumerate(list_of_number_of_edge):
    edges_ = edges[:, :number_of_edge][..., np.newaxis]
    image_rec = mp.dewhitening(mp.reconstruct(edges_, mask=True))
    fig, axs[i_ax] = mp.imshow(image_rec/vmax, fig=fig, ax=axs[i_ax], norm=False, mask=True)
    axs[i_ax].text(5, 29, 'N=%d' % number_of_edge, color='red', fontsize=24)
plt.tight_layout()
fig.subplots_adjust(hspace = .0, wspace = .0, left=0.0, bottom=0., right=1., top=1.)

mp.savefig(fig, name + '_movie')

Running simulations on a set of natural images¶

import pickle with open('mat/EUVIP_sparseness_serre07_distractors_images.pickle', 'rb') as f: imagelist = pickle.load(f) print(imagelist, len(imagelist)) imagelist_ = [] count = 0 for filename, croparea in (imagelist): matname = os.path.join('mat/edges/EUVIP_sparseness_serre07_distractors', filename.replace('.png', '') + str(croparea) + '.npy') if os.path.isfile(matname): count += 1 imagelist_.append([filename, croparea]) if count>=100: break print(len(imagelist_)) pickle.dump(imagelist_, open('mat/EUVIP_sparseness_serre07_distractors_images.pickle', "wb" ) )np.load('mat/edges/' + experiment + '_'+ name_database + '_edges.npy').shape, 2**11

In [16]:

%%writefile experiment_sparseness.py
# -*- coding: utf8 -*-
from __future__ import division, print_function
"""

$ python experiment_sparseness.py

to remove the cached files:
rm -fr **/SparseLets* **/**/SparseLets*

"""
import sys
experiment = sys.argv[1]
parameter_file = sys.argv[2]
name_database = sys.argv[3]
N_image = int(sys.argv[4])
print('N_image', N_image)
N = int(sys.argv[5])
do_linear = (sys.argv[6] == 'True')

import numpy as np
from SparseEdges import SparseEdges
mps = []
for name_database in [name_database]:
    mp = SparseEdges(parameter_file)
    mp.pe.figpath, mp.pe.formats, mp.pe.dpi = 'figures', ['png', 'pdf', 'jpg'], 450
    mp.pe.datapath = 'database/'
    mp.pe.N_image = N_image
    mp.pe.do_mask = True
    mp.pe.N = N
    mp.pe.n_theta = 1
    mp.pe.B_theta = np.inf
    mp.init()
    # normal experiment
    imageslist, edgeslist, RMSE = mp.process(exp=experiment, name_database=name_database)
    mps.append(mp)
    # control experiment
    if do_linear:
        mp.pe.MP_alpha = np.inf
        mp.init()
        imageslist, edgeslist, RMSE = mp.process(exp=experiment + '_linear', name_database=name_database)
        mps.append(mp)

Overwriting experiment_sparseness.py

In [17]:

if cluster:
    for cmd in [
        "frioul_list_jobs -v |grep job_array_id |uniq -c",
    ]:
        print(run_on_cluster(cmd))

In [18]:

experiment_folder = experiment = 'retina-sparseness'

cluster = True
cluster = False

do_update = True
do_update = False

do_cleanup = False
do_cleanup = True

do_run = False
do_run = True

experiments = [experiment]

def run_cmd(cmd, doit=True):
    import subprocess
    print ('⚡︎ Running ⚡︎ ', cmd)
    if doit:
        stdout = subprocess.check_output([cmd], shell=True)
        return stdout.decode()#.splitlines()

SERVER = 'perrinet.l@frioul.int.univ-amu.fr'
PATH = '/hpc/invibe/perrinet.l/science/{}/'.format(experiment_folder)

def push_to_cluster(source="{results,data_cache,experiment_sparseness.py,database}", 
                       PATH=PATH, SERVER=SERVER, 
                       opts="-av -u --exclude .AppleDouble --exclude .git"):
    fullcmd = 'ssh {} "mkdir -p {} " ; '.format(SERVER, PATH)
    fullcmd += 'rsync {} {}  {}:{} '.format(opts, source, SERVER, PATH)
    return run_cmd (fullcmd)

def run_on_cluster(cmd, PATH=PATH, SERVER=SERVER):
    import subprocess
    fullcmd = 'ssh {SERVER} "cd {PATH} ; {cmd} "'.format(SERVER=SERVER, PATH=PATH, cmd=cmd)
    return run_cmd (fullcmd)

def pull_from_cluster(source="{results,data_cache,debug.log}", dest=".", 
                       PATH=PATH, SERVER=SERVER, 
                       opts="-av -u --delete --exclude .AppleDouble --exclude .git"):
    fullcmd = 'rsync {} {}:{}{} {} '.format(opts, SERVER, PATH, source, dest)
    return run_cmd (fullcmd)

# update
if cluster and do_update:
    print(run_on_cluster("frioul_batch  'cd /hpc/invibe/perrinet.l/science/SparseEdges/; make update_dev'"))

# clean-up
if cluster and do_cleanup:
    push_to_cluster()

    for cmd in [
        #"rm -fr results data_cache ",
        "find . -name *lock* -exec rm -fr {} \\;",
        "touch frioul; rm frioul* ",
        ]:
        print(run_on_cluster(cmd))

# RUNNING
if do_run:
    if cluster:
        fullcmd = 'ipython  experiment_sparseness.py {experiment} {parameter_file} {name_database} {N_image} {N} {do_linear} '.format(
            experiment=experiment, parameter_file=parameter_file, 
            name_database=name_database, N_image=N_image, N=N, do_linear=do_linear)
        for cmd in [
            "frioul_batch  -M 136 '{}' ".format(fullcmd), 
            "frioul_list_jobs -v |grep job_array_id |uniq -c",
                    ]:
            print(run_on_cluster(cmd))
    else:
        fullcmd = 'ipython3 experiment_sparseness.py {experiment} {parameter_file} {name_database} {N_image} {N} {do_linear} '.format(
            experiment=experiment, parameter_file=parameter_file, 
            name_database=name_database, N_image=N_image, N=N, do_linear=do_linear)
        run_cmd (fullcmd)

⚡︎ Running ⚡︎  ipython3 experiment_sparseness.py retina-sparseness https://raw.githubusercontent.com/bicv/SparseEdges/master/default_param.py serre07_distractors 100 2048 False

In [19]:

import time, os
# GETTING the data
import time, os
while True:
    if cluster:
        print(pull_from_cluster())
        print(run_on_cluster("tail -n 10 {}".format(os.path.join(PATH, 'debug.log'))))
        print(run_on_cluster("frioul_list_jobs -v |grep job_array_id |uniq -c"))
    locks = run_cmd ("find . -name *lock -exec ls -l {} \;")
    print(locks)
    if len(locks) == 0: break
    time.sleep(100)

⚡︎ Running ⚡︎  find . -name *lock -exec ls -l {} \;

In [20]:

!ssh perrinet.l@frioul.int.univ-amu.fr "python -c'import numpy as np; print(np.pi)'"

3.14159265359

In [21]:

%%bash
ssh perrinet.l@frioul.int.univ-amu.fr "python -c'import numpy as np; print(np.pi)'"

3.14159265359

Analysing results¶

First, we retrieve edges from a prior edge extraction

name_database='serre07_distractors'¶

imageslist, edgeslist, RMSE = mp.process(exp=experiment, name_database=name_database)

In [22]:

%run  experiment_sparseness.py retina-sparseness https://raw.githubusercontent.com/bicv/SparseEdges/master/default_param.py serre07_distractors 100 2048 False

N_image 100

In [23]:

imageslist, edgeslist, RMSE = mp.process(exp=experiment, name_database=name_database)

In [24]:

edgeslist

Out[24]:

array([[[ 140.  ,  122.  ,  140.  , ...,  168.  ,  163.  ,   93.  ],
        [ 141.  ,  113.  ,  148.  , ...,  168.  ,  163.  ,  146.  ],
        [ 140.  ,  194.  ,  151.  , ...,  168.  ,  166.  ,   57.  ],
        ..., 
        [  94.  ,  134.  ,   55.  , ...,  187.  ,  213.  ,  203.  ],
        [ 169.  ,   87.  ,   44.  , ...,  176.  ,  144.  ,  203.  ],
        [ 170.  ,   88.  ,   44.  , ...,  101.  ,  154.  ,  203.  ]],

       [[ 103.  ,   71.  ,  196.  , ...,   85.  ,  118.  ,  151.  ],
        [ 101.  ,  139.  ,  193.  , ...,   82.  ,  118.  ,  102.  ],
        [  98.  ,  151.  ,  194.  , ...,   82.  ,  123.  ,  160.  ],
        ..., 
        [ 167.  ,  177.  ,  178.  , ...,   66.  ,   89.  ,  101.  ],
        [  81.  ,  168.  ,   96.  , ...,   69.  ,   47.  ,  101.  ],
        [  88.  ,  164.  ,   96.  , ...,  107.  ,   48.  ,  102.  ]],

       [[   1.57,    1.57,    1.57, ...,    1.57,    1.57,    1.57],
        [   1.57,    1.57,    1.57, ...,    1.57,    1.57,    1.57],
        [   1.57,    1.57,    1.57, ...,    1.57,    1.57,    1.57],
        ..., 
        [   1.57,    1.57,    1.57, ...,    1.57,    1.57,    1.57],
        [   1.57,    1.57,    1.57, ...,    1.57,    1.57,    1.57],
        [   1.57,    1.57,    1.57, ...,    1.57,    1.57,    1.57]],

       [[   0.24,    0.24,    0.24, ...,    0.24,    0.38,    0.24],
        [   0.24,    0.24,    0.24, ...,    0.24,    0.62,    0.24],
        [   0.15,    0.24,    0.24, ...,    0.38,    0.24,    0.15],
        ..., 
        [   0.24,    0.38,    0.24, ...,    0.38,    0.38,    0.38],
        [   0.02,    0.62,    0.62, ...,    0.62,    0.06,    0.24],
        [   0.09,    0.38,    0.24, ...,    0.24,    0.06,    0.24]],

       [[   1.47,    2.43,    1.67, ...,    2.01,    1.48,    1.13],
        [   1.52,    2.26,    1.66, ...,    2.21,    1.57,    0.99],
        [   1.39,    2.05,    1.61, ...,    2.16,    1.41,    0.93],
        ..., 
        [   0.12,    0.45,    0.34, ...,    0.33,    0.24,    0.2 ],
        [   0.12,    0.43,    0.34, ...,    0.31,    0.24,    0.18],
        [   0.12,    0.45,    0.35, ...,    0.31,    0.25,    0.19]],

       [[  -0.  ,    0.  ,   -0.  , ...,    3.14,   -3.14,   -3.14],
        [  -0.  ,    0.  ,   -0.  , ...,    3.14,    0.  ,   -3.14],
        [  -0.  ,   -3.14,    0.  , ...,   -0.  ,    0.  ,   -0.  ],
        ..., 
        [  -3.14,   -3.14,   -3.14, ...,   -0.  ,   -3.14,   -0.  ],
        [   3.14,    3.14,   -0.  , ...,    3.14,    3.14,   -3.14],
        [  -0.  ,    3.14,    3.14, ...,   -0.  ,    3.14,   -0.  ]]])

In [25]:

fig, [A, B] = plt.subplots(1, 2, figsize=(fig_width, fig_width/1.618), subplot_kw={'axisbg':'w'})
A.set_color_cycle(np.array([[1., 0., 0.]]))
imagelist, edgeslist, RMSE = mp.process(exp=experiment, name_database=name_database)
RMSE /= RMSE[:, 0][:, np.newaxis]
#print( RMSE.shape, edgeslist.shape)
value = edgeslist[4, ...]
#value /= value[0, :][np.newaxis, :]
value /= RMSE[:, 0][np.newaxis, :]

B.semilogx( value, alpha=.7)

A.semilogx( RMSE.T, alpha=.7)
A.set_xlabel('l0')
B.set_xlabel('l0')
A.axis('tight')
B.axis('tight')
_ = A.set_ylabel('RMSE')


#plt.locator_params(axis = 'x', nbins = 5)
#plt.locator_params(axis = 'y', nbins = 5)

mp.savefig(fig, experiment + '_raw')

In [26]:

imagelist, edgeslist, RMSE = mp.process(exp=experiment + '_linear', name_database=name_database)
RMSE /= RMSE[:, 0][:, np.newaxis]
print(RMSE, RMSE.shape, edgeslist.shape)

Checking dependence in serre07_distractors_retina-sparseness_linear
------------------------------------------------------------
Entropy: 0.920346740735
------------------------------------------------------------
------------------------------------------------------------
[['d', 'phi', 'theta', 'scale']] KL= 0.00000 ; 0.000
[['phi', 'theta', 'scale'], ['d']] KL= 0.05398 ; 5.865
[['theta', 'scale', 'd'], ['phi']] KL= 0.01775 ; 1.929
[['scale', 'd', 'phi'], ['theta']] KL= -0.00000 ; -0.000
[['d', 'phi', 'theta'], ['scale']] KL= 0.06835 ; 7.427
[['phi', 'theta'], ['scale', 'd']] KL= 0.01775 ; 1.929
[['theta', 'scale'], ['d', 'phi']] KL= 0.06835 ; 7.427
[['scale', 'phi'], ['d', 'theta']] KL= 0.05398 ; 5.865
[['phi', 'theta'], ['scale'], ['d']] KL= 0.07015 ; 7.622
[['theta', 'scale'], ['d'], ['phi']] KL= 0.07015 ; 7.622
[['scale', 'd'], ['phi'], ['theta']] KL= 0.01775 ; 1.929
[['d', 'phi'], ['theta'], ['scale']] KL= 0.06835 ; 7.427
[['d'], ['phi'], ['theta'], ['scale']] KL= 0.07015 ; 7.622
------------------------------------------------------------

[[ nan  nan  nan ...,  nan  nan  nan]
 [ nan  nan  nan ...,  nan  nan  nan]
 [ nan  nan  nan ...,  nan  nan  nan]
 ..., 
 [ nan  nan  nan ...,  nan  nan  nan]
 [ nan  nan  nan ...,  nan  nan  nan]
 [ nan  nan  nan ...,  nan  nan  nan]] (100, 4096) (6, 4096, 100)

In [27]:

fig = plt.figure(figsize=(fig_width/1.618, fig_width/1.618))
if do_linear:
    fig, ax, inset = mp.plot(mps=[mp, mp], experiments=[experiment, experiment + '_linear'], 
                             databases=[name_database, name_database], fig=fig, 
                             color=[0., 0., 1.], scale=False, labels=['MP', 'lin'])
else:
    fig, ax, inset = mp.plot(mps=[mp], experiments=[experiment], databases=[name_database], fig=fig, 
                  color=[0., 0., 1.], scale=False, labels=['MP'])

mp.savefig(fig, experiment + '_raw_inset')

trying different fits¶

on the modulation function¶

!pip install lmfit

from lmfit.models import ExpressionModel import numpy as np mod = ExpressionModel('off + amp * exp(-x/x0) * sin(x*phase)') x = np.linspace(0, 10, 501) params = mod.make_params(off=0.25, amp=1.0, x0=2.0, phase=0.04) y = mod.eval(params, x=x) yfig, A = plt.subplots(1, 1, figsize=(fig_width, fig_width/1.618), subplot_kw={'axisbg':'w'}) from lmfit.models import ExpressionModel #mod = PowerLawModel() mod = ExpressionModel('rho**x') RMSE /= RMSE[:, 0][:, np.newaxis] N = RMSE.shape[1] #number of edges rho = np.zeros(RMSE.shape[0]) for i_image in range(RMSE.shape[0]): #pars = mod.guess(RMSE[i_image, :], x=np.arange(N)) mod.def_vals = {'rho':.99} out = mod.fit(RMSE[i_image, :], x=np.arange(N), verbose=False, weights=1/(np.arange(N)+1)) #print(out.fit_report(min_correl=0.25)) rho[i_image] = out.params.get('rho').value #print 'rho=', rho[i_image] #N_theta = np.log(threshold)/np.log(rho) #print N_theta A.semilogx( RMSE[i_image, :], alpha=.7) params = mod.make_params(rho=rho[i_image]) A.semilogx(mod.eval(params, x=np.arange(N)), 'r--', alpha=.7) A.set_xlabel('l0') A.axis('tight') A.axis('tight') _ = A.set_ylabel('RMSE') print ('rho=', rho.mean(), ', +/- ', rho.std())fig, A = plt.subplots(1, 1, figsize=(fig_width, fig_width/1.618), subplot_kw={'axisbg':'w'}) from lmfit.models import ExpressionModel #mod = PowerLawModel() mod = ExpressionModel('1 - (1- eps_inf) * ( 1 - rho**(x+1))') threshold = .5 verbose = False imagelist, edgeslist, RMSE = mp.process(exp=experiment, name_database=name_database) RMSE /= RMSE[:, 0][:, np.newaxis] N = RMSE.shape[1] #number of edges for i_image in range(RMSE.shape[0]): #pars = mod.guess(RMSE[i_image, :], x=np.arange(N)) mod.def_vals = {'eps_inf':.001, 'rho':.99} out = mod.fit(RMSE[i_image, :], x=np.arange(N), verbose=verbose, weights=1/(np.arange(N)+1)) #print(out.fit_report(min_correl=0.25)) eps_inf = out.params.get('eps_inf').value rho = out.params.get('rho').value #print rho, eps_inf N_theta = np.log((threshold-eps_inf)/(1-eps_inf))/np.log(rho) #print N_theta A.semilogx( RMSE[i_image, :], alpha=.7) params = mod.make_params(eps_inf=eps_inf, rho=rho) A.semilogx(mod.eval(params, x=np.arange(N)), 'r--', alpha=.7) A.set_xlabel('l0') A.axis('tight') A.axis('tight') _ = A.set_ylabel('RMSE') imagelist, edgeslist, RMSE = mp.process(exp=experiment, name_database=name_database) value = edgeslist[4, ...] value /= RMSE[:, 0][np.newaxis, :] #value /= RMSE[:, 0][np.newaxis, :] #RMSE /= RMSE[:, 0][:, np.newaxis] N_image, N = RMSE.shape #number of images x edges value = value.T from lmfit.models import ExpressionModel mod = ExpressionModel('amplitude * exp ( - .5 * log(x+off)**2 / log(rho+1) **2 )') verbose = False amplitude, rho = np.zeros(N_image), np.zeros(N_image) for i_image in range(N_image): #pars = mod.guess(RMSE[i_image, :], x=np.arange(N)) mod.def_vals = {'amplitude':.01, 'rho':100, 'off':1} params = mod.make_params() out = mod.fit(value[i_image, :], x=np.arange(N), verbose=verbose) #print(out.fit_report()) amplitude[i_image] = out.params.get('amplitude').value rho[i_image] = out.params.get('rho').value off[i_image] = out.params.get('off').value print(off)

In [28]:

imagelist, edgeslist, RMSE = mp.process(exp=experiment, name_database=name_database)
value = edgeslist[4, ...].T
#value /= RMSE[:, 0][np.newaxis, :]
value /= RMSE[:, 0][:, np.newaxis]
#RMSE /= RMSE[:, 0][:, np.newaxis]
N_image, N = RMSE.shape #number of images x edges
#value = value.T

from lmfit.models import ExpressionModel mod = ExpressionModel('amplitude * exp ( - .5 * log(x+1)**2 / log(rho+1) **2 )') verbose = False amplitude, rho = np.zeros(N_image), np.zeros(N_image) for i_image in range(N_image): #pars = mod.guess(RMSE[i_image, :], x=np.arange(N)) mod.def_vals = {'amplitude':.01, 'rho':100} params = mod.make_params() out = mod.fit(value[i_image, :], x=np.arange(N), verbose=verbose) #print(out.fit_report()) amplitude[i_image] = out.params.get('amplitude').value rho[i_image] = out.params.get('rho').valueamplitude, rho

In [29]:

imagelist, edgeslist, RMSE = mp.process(exp=experiment, name_database=name_database)
value = edgeslist[4, ...]
value /= RMSE[:, 0][np.newaxis, :]
#RMSE /= RMSE[:, 0][:, np.newaxis]
N = RMSE.shape[1] #number of edges
value = value.T
print(value.shape, RMSE.shape)

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/1.618), subplot_kw={'axisbg':'w'})
from lmfit.models import ExpressionModel
mod = ExpressionModel('amplitude * exp ( - .5 * log(x+1)**2 / rho **2 )')
verbose = False
amplitude, rho = np.zeros(N_image), np.zeros(N_image)

for i_image in range(RMSE.shape[0]):
    #pars = mod.guess(RMSE[i_image, :], x=np.arange(N))
    mod.def_vals = {'amplitude':.01, 'rho':100}
    params = mod.make_params()
    out  = mod.fit(value[i_image, :], x=np.arange(N), verbose=verbose, params=params)#, weights=np.exp(- np.arange(N) / 200))
    #print(out.params)
    #print(out.fit_report())
    amplitude[i_image] = out.params.get('amplitude').value
    rho[i_image] =  out.params.get('rho').value

    ax.loglog( value[i_image, :], alpha=.2)
    params = mod.make_params(amplitude=amplitude[i_image], rho=rho[i_image])
    ax.loglog(mod.eval(params, x=np.arange(N)), 'r--', alpha=.2)
    ax.set_xlabel('l0')
    ax.axis('tight')
    _ = ax.set_ylabel('coefficient')

(100, 4096) (100, 4096)

In [30]:

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/1.618), subplot_kw={'axisbg':'w'})

for i_image in range(N_image):
    ax.loglog( value[i_image, :], alpha=.2)
    params = mod.make_params(amplitude=amplitude[i_image], rho=rho[i_image])
    ax.loglog(mod.eval(params, x=np.arange(N)), 'r--', alpha=.2)
    ax.set_xlabel('l0')
    ax.axis('tight')
    _ = ax.set_ylabel('coefficient')            
mp.savefig(fig, experiment + '_fit_all')

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/1.618), subplot_kw={'axisbg':'w'}) for i_image in range(N_image): ax.semilogy(np.log((np.arange(N)+1)/np.log(rho[i_image]+1)), value[i_image, :]/amplitude[i_image], alpha=.2) params = mod.make_params(amplitude=amplitude[i_image], rho=rho[i_image]) ax.semilogy(np.log((np.arange(N)+1)/np.log(rho[i_image]+1)), mod.eval(params, x=np.arange(N))/amplitude[i_image], 'r--', alpha=.2) ax.set_xlabel('l0 norm') ax.axis('tight') _ = ax.set_ylabel('norm. coefficient') mp.savefig(fig, experiment + '_fit_norm')

In [31]:

fig, axs = plt.subplots(1, 3, figsize=(fig_width, fig_width/1.618), subplot_kw={'axisbg':'w'})

axs[0].hist(amplitude)
axs[1].hist(np.abs(rho))
axs[2].scatter(amplitude, np.abs(rho))
for ax in axs: 
    ax.axis('tight')
    _ = ax.set_ylabel('')            
    _ = ax.set_yticks([])            
axs[0].set_ylabel('probability')            
axs[0].set_xlabel('amplitude')
axs[1].set_xlabel('rho')
axs[2].set_xlabel('amplitude')
axs[2].set_ylabel('rho')
fig.tight_layout()
mp.savefig(fig, experiment + '_fit_hist')

In [33]:

fig, axs = plt.subplots(1, 1, figsize=(fig_width/2.618, fig_width/1.618), subplot_kw={'axisbg':'w'})

axs.hist(np.abs(rho))
axs.axis('tight')
_ = axs.set_ylabel('')            
_ = axs.set_yticks([])            
axs.set_ylabel('probability')            
axs.set_xlabel('amplitude')
fig.tight_layout()
mp.savefig(fig, experiment + '_fit_hist')

[autoreload of imageio.plugins.pillow failed: Traceback (most recent call last):
  File "/usr/local/lib/python3.6/site-packages/IPython/extensions/autoreload.py", line 247, in check
    superreload(m, reload, self.old_objects)
ValueError: A Format named 'BMP-PIL' is already registered, use overwrite=True to replace.
]

on the pdf¶

help(np.histogram)

In [34]:

value.max(axis=1).shape

Out[34]:

(100,)

In [35]:

%pwd

Out[35]:

'/Users/lolo/pool/science/RetinaClouds/2016-12-25_DropLets_ms'

In [36]:

#imagelist, edgeslist, RMSE = mp.process(exp=experiment + '_linear', name_database=name_database)
#imagelist, edgeslist, RMSE = mp.process(exp=experiment, name_database=name_database)
edgeslist = np.load('data_cache/edges/' + experiment + '_' + name_database + '_edges.npy')
value = edgeslist[4, ...].T
#value /= RMSE[:, 0][np.newaxis, :]
value /= value.min(axis=1)[:, np.newaxis]
#RMSE /= RMSE[:, 0][:, np.newaxis]
N_image, N = value.shape #number of images x edges
#value = value.T

In [37]:

N_bins, a_max = 128, value.max()
start, end = N_bins//16, N_bins
print(a_max)
v_hist = np.zeros((N_image, N_bins))
#bins = np.linspace(0, a_max, N_bins+1, endpoint=True)#[:-1]
#print(bins.shape)
for i_image in range(N_image):
    #v_hist[i_image, : ], v_bins = np.histogram(value[i_image, :], bins=bins) 
    v_hist[i_image, : ], v_bins = np.histogram(value[i_image, :], bins=N_bins) 
    v_hist[i_image, : ] /= v_hist[i_image, : ].sum()
print(v_bins.shape)
v_middle = .5*(v_bins[1:]+v_bins[:-1])

12.9868477695
(129,)

In [38]:

plt.plot(v_bins[1:], v_middle)
print(v_bins[0], v_middle[0])

print(v_bins[0], v_middle[0])

print(start, end)

1.0 1.02124006989
1.0 1.02124006989
8 128

MLE estimate of rho: https://en.wikipedia.org/wiki/Power_law#Maximum_likelihood

In [39]:

amplitude, rho = np.zeros(N_image), np.zeros(N_image)
for i_image in range(N_image):
    rho[i_image] =  1 +  (end-start) / np.sum(np.log(value[i_image, start:end]))
    amplitude[i_image] = rho[i_image] - 1
print(rho)

[ 1.48  1.77  1.86  1.75  1.67  1.86  1.85  1.6   1.77  2.15  1.73  1.92
  1.73  1.63  1.72  1.77  1.71  1.73  1.71  1.52  1.86  1.84  2.01  1.9
  1.67  1.69  2.07  1.73  1.61  1.7   1.56  1.79  1.69  1.74  1.81  1.7
  1.66  1.82  1.98  1.7   1.75  1.6   1.86  1.7   1.73  1.81  1.94  1.65
  1.67  1.6   1.56  1.72  1.84  1.81  1.88  1.93  1.95  1.69  1.65  1.66
  1.69  1.53  1.67  1.72  1.75  1.83  1.72  1.73  1.53  1.61  1.63  1.52
  1.83  1.73  1.66  1.76  1.84  1.77  1.75  1.82  1.71  1.99  1.76  1.73
  1.72  1.73  1.92  1.95  1.53  1.85  1.73  1.9   1.77  1.59  1.74  1.68
  1.93  1.64  1.7   1.74]

In [40]:

import lmfit
from lmfit import Model

def model(x, A, x_0, B):
    f =  A / x * np.exp( -.5 * np.log(x/x_0)**2 / B**2 )
    #f /= f.sum()
    return f

weights = np.linspace(0, 1, N_bins)
weights = np.linspace(1, 0, N_bins)
weights = np.ones(N_bins)
verbose = False

A, x_0, B = np.zeros(N_image), np.zeros(N_image), np.zeros(N_image)
for i_image in range(N_image):
    mod = Model(model)
    mod.set_param_hint('A', value=.05, min=0.)
    #mod.set_param_hint('x_0', value=.45, min=0.45, max=0.46)
    mod.set_param_hint('x_0', value=.5, min=0.)
    mod.set_param_hint('B', value=1.9, min=0.)

    valid = (v_hist[i_image, :] > 0.)
    out  = mod.fit(v_hist[i_image, valid], x=v_middle[valid], 
                   verbose=verbose, weights=weights[valid], method='leastsq', maxfev=1e6)
    if verbose: print(out.fit_report())
    A[i_image] = out.params.get('A').value
    x_0[i_image] =  out.params.get('x_0').value
    B[i_image] =  out.params.get('B').value
print ('A=', A.mean(), ', +/- ', A.std())    
print ('x_0=', x_0.mean(), ', +/- ', x_0.std())    
print ('B=', B.mean(), ', +/- ', B.std())    

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/1.618), subplot_kw={'axisbg':'w'})
for i_image in range(N_image):
    ax.plot(v_middle[valid], v_hist[i_image, valid], '.', alpha=.2)
    valid = (v_hist[i_image, :] > 0.)
    #params = mod.make_params(A=A[i_image], x_0=x_0[i_image], B=B[i_image])
    #ax.plot(v_middle[valid], mod.eval(params, x=v_middle[valid]), 'r', alpha=.2)
    ax.plot(v_middle[valid], model(v_middle[valid], A=A[i_image], x_0=x_0[i_image], B=B[i_image]), 'r', alpha=.2)
    ax.set_yscale('log')
    ax.set_xscale('log')
    ax.axis('tight')
    #ax.set_xlim(a_min, a_max)
    ax.set_ylim(.0003, .1)
    ax.set_ylabel('density')
    ax.set_xlabel('coefficient')            
mp.savefig(fig, experiment + '_proba')

A= 0.0576047499931 , +/-  0.01380543006
x_0= 1.19630015509 , +/-  0.178766832634
B= 0.475070357045 , +/-  0.063421389447

In [41]:

import lmfit
from lmfit import Model

def model(x, A, rho):
    f =  A / x ** rho
    #f /= f.sum()
    return f

weights = np.linspace(1, 0, N_bins)
weights = np.linspace(0, 1, N_bins)
weights = np.ones(N_bins)
verbose = False

A, rho = np.zeros(N_image), np.zeros(N_image)
for i_image in range(N_image):
    mod = Model(model)
    mod.set_param_hint('A', value=.05, min=0.)
    mod.set_param_hint('rho', value=2.5, min=1.)

    valid = (v_hist[i_image, :] > 0.)
    out  = mod.fit(v_hist[i_image, valid], x=v_middle[valid], 
                   verbose=verbose, weights=weights[valid], method='leastsq', maxfev=1e6)
    if verbose: print(out.fit_report())
    A[i_image] = out.params.get('A').value
    rho[i_image] =  out.params.get('rho').value
print ('A=', A.mean(), ', +/- ', A.std())    
print ('rho=', rho.mean(), ', +/- ', rho.std())    

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/1.618), subplot_kw={'axisbg':'w'})
for i_image in range(N_image):
    ax.plot(v_middle[valid], v_hist[i_image, valid], '.', alpha=.2)
    valid = (v_hist[i_image, :] > 0.)
    ax.plot(v_middle[valid], model(v_middle[valid], A=A[i_image], rho=rho[i_image]), 'r', alpha=.2)
    ax.set_yscale('log')
    ax.set_xscale('log')
    ax.axis('tight')
    ax.set_ylabel('density')
    ax.set_xlabel('coefficient')            
mp.savefig(fig, experiment + '_proba')

A= 0.0662767889623 , +/-  0.0156305065431
rho= 2.24177826437 , +/-  0.406469456929

In [42]:

from lmfit.models import ExpressionModel
mod = ExpressionModel('amplitude * x**-rho ')
#mod = ExpressionModel('amplitude * exp( - log(x)**2/rho**2 ) ')
#mod = ExpressionModel('amplitude * exp( - x/rho ) ')
verbose = False
for i_image in range(N_image):
    #pars = mod.guess(RMSE[i_image, :], x=np.arange(N))
    mod.def_vals = {'amplitude': amplitude[i_image], 'rho': rho[i_image]}
    params = mod.make_params()
    out  = mod.fit(v_hist[i_image, start:end], x=v_middle[start:end], verbose=verbose)
    #print(out.fit_report())
    amplitude[i_image] = out.params.get('amplitude').value
    rho[i_image] =  out.params.get('rho').value
print(rho)

[ 3.2   3.    2.68  2.52  2.89  2.98  2.32  3.31  3.37  2.49  2.85  2.85
  3.01  3.02  4.29  2.68  2.83  2.92  2.86  3.37  2.41  2.7   2.3   3.2
  2.88  3.86  2.31  2.93  3.27  4.77  2.51  2.85  2.83  2.38  2.46  3.28
  3.11  3.05  2.93  2.66  2.79  2.97  2.9   2.52  2.47  3.8   2.67  3.19
  3.16  2.54  2.97  3.11  2.26  2.88  2.38  2.3   2.89  3.27  2.97  2.73
  3.52  3.74  2.62  2.7   2.77  2.8   3.41  3.02  3.17  3.7   3.43  4.05
  2.82  3.29  2.96  3.03  2.76  3.47  2.96  3.6   3.56  2.13  2.9   2.72
  2.62  3.21  3.16  2.69  3.83  2.9   2.44  2.66  3.31  3.7   3.12  2.94
  2.39  3.07  3.07  3.02]

In [43]:

fig, ax = plt.subplots(1, 1, figsize=(fig_width/3, fig_width/3), subplot_kw={'axisbg':'w'})

for i_image in range(N_image):
    ax.plot(v_middle, v_hist[i_image, :], alpha=.2)
    params = mod.make_params(amplitude=amplitude[i_image], rho=rho[i_image])
    ax.plot(v_middle[start:end], mod.eval(params, x=v_middle[start:end]), 'r.', alpha=.2)
    if True:
      ax.set_yscale('log')
      ax.set_xscale('log')
    ax.axis('tight')
    ax.set_xlim(1.5, 5)
    ax.set_ylim(.0003, .05)
    ax.set_ylabel('density')
    ax.set_xlabel('coefficient')            
mp.savefig(fig, experiment + '_proba')

In [44]:

fig, axs = plt.subplots(1, 1, figsize=(fig_width/3, fig_width/3), subplot_kw={'axisbg':'w'})

axs.hist(np.abs(rho), bins=np.linspace(2, 4, 5))
axs.axis('tight')
_ = axs.set_ylabel('')            
_ = axs.set_yticks([])            
axs.set_ylabel('probability density')            
axs.set_xlabel(r'$\rho$')
fig.tight_layout()
mp.savefig(fig, experiment + '_fit_hist')

In [45]:

#fig, ax = plt.subplots(1, 1, figsize=(fig_width/2, fig_width/2), subplot_kw={'axisbg':'w'})
fig = plt.figure(figsize=(fig_width/2, fig_width/2))

ax = fig.add_axes([0.18, 0.15, .8, .8], axisbg='w')

for i_image in range(N_image):
    ax.plot(v_middle, v_hist[i_image, :], '-', alpha=.1, lw=.5)
    params = mod.make_params(amplitude=amplitude[i_image], rho=rho[i_image])
    ax.plot(v_middle[start:end], mod.eval(params, x=v_middle[start:end]), 'r', alpha=.2, lw=.5)
    if True:
      ax.set_yscale('log')
      ax.set_xscale('log')
    ax.axis('tight')
    ax.set_xlim(1.5, 5.)
    ax.set_ylim(.003, .05)
    ax.set_ylabel('probability density')
    ax.set_xlabel('coefficient')  

inset = fig.add_axes([0.58, 0.55, .4, .4], axisbg='w')

inset.hist(np.abs(rho))
inset.axis('tight')
_ = inset.set_ylabel('')            
_ = inset.set_yticks([])            
inset.set_ylabel('# occurences')            
inset.set_xlabel(r'$\rho$')

#fig.subplots_adjust(left=0.22, bottom=0.1, right=.9, top=.9)

mp.savefig(fig, experiment + '_proba_inset')

TODO : do the same for experiment + '_linear'¶

In [46]:

rho_0 = rho.mean()
print(rho_0)
v_hist_scale = np.zeros((N_image, N_bins))
for i_image in range(N_image):
    #v_hist[i_image, : ], v_bins = np.histogram(value[i_image, :], bins=bins) 
    v_hist_scale[i_image, : ], v_bins = np.histogram(value[i_image, :]**((rho_0-1)/(rho[i_image]-1)), bins=N_bins) 
    v_hist_scale[i_image, : ] /= v_hist_scale[i_image, : ].sum()

2.98177377922

In [47]:

amplitude_scale, rho_scale = np.zeros(N_image), np.zeros(N_image)
for i_image in range(N_image):
    mod.def_vals = {'amplitude': amplitude[i_image], 'rho': rho[i_image]}
    params = mod.make_params()
    out  = mod.fit(v_hist_scale[i_image, start:end], x=v_middle[start:end], verbose=verbose)
    amplitude_scale[i_image] = out.params.get('amplitude').value
    rho_scale[i_image] =  out.params.get('rho').value
print(rho_scale)

[ 3.    2.99  3.02  2.99  2.97  2.99  2.96  3.03  3.2   2.91  2.97  2.95
  2.97  2.99  3.08  2.99  3.03  2.95  3.07  3.26  2.99  3.    2.92  3.09
  2.92  3.35  2.89  2.96  3.05  3.09  2.95  2.96  2.95  3.    3.    3.01
  2.94  2.98  3.    2.95  2.97  2.99  2.97  2.91  2.94  2.94  2.89  2.92
  3.05  2.88  2.99  3.03  2.92  2.91  2.93  2.9   2.98  3.03  2.96  2.89
  3.03  3.15  2.92  2.89  2.89  3.    3.06  2.98  2.95  2.97  3.07  3.15
  2.94  3.03  2.96  2.96  2.97  2.98  2.97  3.13  3.02  2.93  3.01  3.03
  2.95  3.05  3.01  2.86  2.96  3.01  2.86  2.88  3.02  3.16  2.96  3.01
  2.95  3.    3.05  3.02]

In [48]:

fig, ax = plt.subplots(1, 1, figsize=(fig_width, fig_width/1.618), subplot_kw={'axisbg':'w'})

for i_image in range(N_image):
    ax.plot(v_middle, v_hist_scale[i_image, :], alpha=.2)
    params = mod.make_params(amplitude=amplitude_scale[i_image], rho=rho_scale[i_image])
    ax.plot(v_middle[start:end], mod.eval(params, x=v_middle[start:end]), 'r.', alpha=.2)
    if True:
      ax.set_yscale('log')
      ax.set_xscale('log')
    ax.set_xlim(1.25, 4)
    ax.set_ylim(.0009, .05)
    ax.set_ylabel('density')
    ax.set_xlabel('coefficient')            
mp.savefig(fig,  experiment + '_proba_scaled')

In [49]:

plt.hist(rho_scale)

Out[49]:

(array([ 10.,  26.,  34.,  17.,   6.,   2.,   2.,   1.,   1.,   1.]),
 array([ 2.86,  2.91,  2.96,  3.01,  3.06,  3.11,  3.15,  3.2 ,  3.25,
         3.3 ,  3.35]),
 <a list of 10 Patch objects>)

DropLets¶

In [50]:

from scipy.stats import powerlaw
N_sparse = 6
sparseness =  np.linspace(2, 7, N_sparse, endpoint=True)
N_edge = N

fig , ax = plt.subplots()
bins= np.linspace(1, 10, 100)
for a in sparseness:
    #s =  np.random.power(a=a, size=N_edge)
    s = 1/powerlaw.rvs(a=a, size = N_edge)
    hist, bins_ = np.histogram(s, bins=bins)
    ax.loglog(bins[1:], hist, label=a)
ax.legend()

Out[50]:

<matplotlib.legend.Legend at 0x10a17fef0>

help(powerlaw )from scipy.stats import powerlaw r = powerlaw.rvs(a = 5, size = 1000) r

In [51]:

frames =[]
for a in sparseness:
    frames.append(mp.texture( N_edge=N_edge, a=a, randn=False))

In [52]:

fig, axs = plt.subplots(1, N_sparse, sharex=True, sharey=True)
fig.set_size_inches(fig_width, 1.2*fig_width/N_sparse)
for i_sparse in range( N_sparse):
    vmax=np.abs(frames[i_sparse]).max()
    vmin=-vmax
    axs[i_sparse].imshow(frames[i_sparse], origin='lower', cmap='gray', vmin=vmin, vmax=vmax, interpolation='none')
    axs[i_sparse].axis('tight')
    axs[i_sparse].set_xticks([])
    axs[i_sparse].set_yticks([])
    axs[i_sparse].set_title(label = r'$\rho=%.0f$' % sparseness[i_sparse])

fig.tight_layout()
fig.subplots_adjust(left=0, bottom=0, right=1, top=.8, wspace=0, hspace=0)   

mp.savefig(fig, 'droplets')

In [5]:

import numpy as np
import MotionClouds as mc
import matplotlib.pyplot as plt

# PARAMETERS
seed = 2042
np.random.seed(seed=seed)
N_sparse = 5
sparse_base = 2.e5
sparseness =  np.logspace(-1, 0, N_sparse, base=sparse_base, endpoint=True)
print(sparseness)

# TEXTON
N_X, N_Y, N_frame = 256, 256, 1
fx, fy, ft = mc.get_grids(N_X, N_Y, 1)
mc_i = mc.envelope_gabor(fx, fy, ft, sf_0=0.05, B_sf=0.025, B_theta=np.inf)
#print(ft.shape)
#print(mc_i.shape)
#fig, axs = plt.subplots(1, 1, figsize=(fig_width, fig_width))
#axs.imshow(mc.envelope_speed(fx, fy, ft)[:, :, 0], vmin=-1, vmax=1, cmap=plt.gray())

#texton = 2*mc.rectif(mc.random_cloud(mc_i, impulse=True))-1
#fig, axs = plt.subplots(1, 1, figsize=(fig_width, fig_width))
#axs.imshow(texton[:, :, 0], vmin=-1, vmax=1, cmap=plt.gray())

values = np.random.randn(N_X, N_Y, N_frame)
#a = 2.
#values = np.random.pareto(a=a, size=(N_X, N_Y, N_frame)) + 1
#values *= np.sign(np.random.randn(N_X, N_Y, N_frame))

#chance = np.random.rand(N_X, N_Y, N_frame)
chance = np.argsort(-np.abs(values.ravel()))
#fig, axs = plt.subplots(1, 1, figsize=(fig_width, fig_width))
#axs.plot(np.abs(values.ravel())[chance])
chance = np.array(chance, dtype=np.float)
chance /= chance.max()
chance = chance.reshape((N_X, N_Y, N_frame))
#print(chance.min(), chance.max())
#fig, axs = plt.subplots(1, 1, figsize=(fig_width, fig_width))
#axs.imshow(chance[:, :, 0], vmin=0, vmax=1, cmap=plt.gray())

fig, axs = plt.subplots(1, N_sparse, figsize=(fig_width, fig_width/N_sparse))
for i_ax, l0_norm in enumerate(sparseness):

    threshold = 1 - l0_norm
    mask = np.zeros_like(chance)
    mask[chance > threshold] = 1.
    
    im = 2*mc.rectif(mc.random_cloud(mc_i, events=mask*values))-1
                
    axs[i_ax].imshow(im[:, :, 0], vmin=-1, vmax=1, cmap=plt.gray())
    #axs[i_ax].text(9, 80, r'$n=%.0f\%%$' % (noise*100), color='white', fontsize=10)
    axs[i_ax].text(4, 40, r'$\epsilon=%.0e$' % l0_norm, color='white', fontsize=8)
    axs[i_ax].set_xticks([])
    axs[i_ax].set_yticks([])
plt.tight_layout()
fig.subplots_adjust(hspace = .0, wspace = .0, left=0.0, bottom=0., right=1., top=1.)
# plt.savefig(fig, experiment + '_droplets')

[5.00000000e-06 1.05737126e-04 2.23606798e-03 4.72870805e-02
 1.00000000e+00]

/var/folders/3p/m0g52j9j69z3gj8ktpgg1dm00000gn/T/ipykernel_67799/2693606862.py:35: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  chance = np.array(chance, dtype=np.float)

modeling different sparseness exponents in a natural image¶

In [ ]:

!ls data_cache/retina-lena.npy

In [ ]:

print(mp.pe.matpath)

In [ ]:

name = experiment.replace('sparseness', 'lena')
matname = os.path.join(mp.pe.matpath, name + '.npy')
N_rho = 3
fig, axs = plt.subplots(1, N_rho, figsize=(fig_width, fig_width/N_rho))
vmax = 1.
for i_ax, rho in enumerate(np.logspace(-1, 1, N_rho, base=2)):
    edges = np.load(matname)
    edges[4, :] = edges[4, :] ** rho
    image_rec = mp.dewhitening(mp.reconstruct(edges, mask=True))      
    fig, axs[i_ax] = mp.imshow(image_rec/vmax, fig=fig, ax=axs[i_ax], norm=False, mask=True)
    axs[i_ax].text(5, 29, r'$\rho=%.1f$' % rho, color='red', fontsize=16)
plt.tight_layout()
fig.subplots_adjust(hspace = .0, wspace = .0, left=0.0, bottom=0., right=1., top=1.)

mp.savefig(fig, name + '_rescale')

some book keeping for the notebook¶

In [12]:

%reload_ext watermark

In [10]:

%watermark -i -h -m -v -p MotionClouds,numpy,SLIP,LogGabor,SparseEdges,matplotlib,scipy,pillow,imageio

Python implementation: CPython
Python version       : 3.9.8
IPython version      : 7.28.0

MotionClouds: 20200212
numpy       : 1.21.4
SLIP        : not installed
LogGabor    : not installed
SparseEdges : not installed
matplotlib  : 3.4.3
scipy       : 1.6.0
pillow      : not installed
imageio     : 2.9.0

Compiler    : Clang 13.0.0 (clang-1300.0.29.3)
OS          : Darwin
Release     : 21.1.0
Machine     : x86_64
Processor   : i386
CPU cores   : 4
Architecture: 64bit

Hostname: ekla

MEUL with a non-parametric homeostasis

Laurent Perrinet

2017-11-07 16:13

In this notebook, we will study how homeostasis (cooperation) may be an essential ingredient to this algorithm working on a winner-take-all basis (competition). This extension has been published as Perrinet, Neural Computation (2010) (see https://laurentperrinet.github.io/publication/perrinet-10-shl ). Compared to other posts, such as this previous post, we improve the code to not depend on any parameter (namely the Cparameter of the rescaling function). For that, we will use a non-parametric approach based on the use of cumulative histograms.

This is joint work with Victor Boutin and Angelo Francisioni. See also the other posts on unsupervised learning.

Designing a A0 poster using matplotlib

Laurent Perrinet

2017-10-25 14:05

Laurent Perrinet

2017-06-15 16:42

Lors de la visite au laboratoire d'une brillante élève de seconde (salut Lena!), nous avons inventé ensemble un jeu: le jeu de l'urne. Le principe est simple: il faut deviner la couleur de la balle qu'on tire d'une urne contenant autant de balles rouges que noires - et ceci le plus tôt possible. Plus précisément, les règles sont:

On a un ensemble de balles, la motié sont rouges, l'autre moitié noires (c'est donc un nombre pair de balles qu'on appelera $N$, disons $N=8$).
Elles sont dans une urne opaque et donc on ne peut pas les voir à moins de les tirer une par une (sans remise dans l'urne). On peut tirer autant de balles qu'on veut pour les observer.
Le but est de deviner la balle qu'on va tirer. Si on gagne (on a bien prédit la couleur), alors on gagne autant de points que le nombre de balles qui étaient dans l'urne au moment de la décision. Sinon on perd autant de points que l'on en aurait gagné!
à long terme, la stratégie du jeu est de décider le meilleur moment où on est prêt à deviner la couleur de la balle qu'on va prendre et ainsi de gagner le plus de points possibles.

Nous avons d'abord créé ce jeu grâce au language de programmation Scratch sur https://scratch.mit.edu/projects/165806365/:

Ici, nous allons essayer de l'analyser plus finement.

testing COMPs-fastPcum_scripted

Laurent Perrinet

2017-05-23 16:13

In this notebook, we will study how homeostasis (cooperation) may be an essential ingredient to this algorithm working on a winner-take-all basis (competition). This extension has been published as Perrinet, Neural Computation (2010) (see https://laurentperrinet.github.io/publication/perrinet-10-shl ). Compared to the previous post, we integrated the faster code to https://github.com/bicv/SHL_scripts.

See also the other posts on unsupervised learning,

This is joint work with Victor Boutin.

testing COMPs-fastPcum

Laurent Perrinet

2017-05-23 16:13

In this notebook, we will study how homeostasis (cooperation) may be an essential ingredient to this algorithm working on a winner-take-all basis (competition). This extension has been published as Perrinet, Neural Computation (2010) (see https://laurentperrinet.github.io/publication/perrinet-10-shl ). Compared to the previous post, we optimize the code to be faster.

See also the other posts on unsupervised learning,

This is joint work with Victor Boutin.

testing COMPs-Pcum

Laurent Perrinet

2017-05-23 15:58

In this notebook, we will study how homeostasis (cooperation) may be an essential ingredient to this algorithm working on a winner-take-all basis (competition). This extension has been published as Perrinet, Neural Computation (2010) (see https://laurentperrinet.github.io/publication/perrinet-10-shl ). In particular, we will show how one can build the non-linear functions based on the activity of each filter and which implement homeostasis.

See also the other posts on unsupervised learning,

This is joint work with Victor Boutin.

A change in the definition of spatial frequency bandwidth?

Laurent Perrinet

2017-05-23 11:43

Since the beginning, we have used a definition of bandwidth in the spatial frequency domain which was quite standard (see supp material for instance):

$$ \mathcal{E}(f; sf_0, B_{sf}) \propto \frac {1}{f} \cdot \exp \left( -.5 \frac {\log( \frac{f}{sf_0} ) ^2} {\log( 1 + \frac {B_sf}{sf_0} )^2 } \right) $$

This is implemented in the folowing code which reads:

env = 1./f_radius*np.exp(-.5*(np.log(f_radius/sf_0)**2)/(np.log((sf_0+B_sf)/sf_0)**2))

However the one implemented in the code looks different (thanks to Kiana for spotting this!), so that one can think that the code is using:

$$ \mathcal{E}(f; sf_0, B_{sf}) \propto \frac {1}{f} \cdot \exp \left( -.5 \frac {\log( \frac{f}{sf_0} ) ^2} {\log(( 1 + \frac {B_sf}{sf_0} )^2 ) } \right) $$

The difference is minimal, yet very important for a correct definition of the bandwidth!

Manipulating speed in motion clouds

Laurent Perrinet

2017-04-10 13:57

An essential dimension of motion is speed. However, this notion is prone to confusions as the speed that has to be measured can be relative to different object. Is it the speed of pixels? The speed of visual objects? We try to distinguish the latter two in this post.

In [4]:

from IPython.display import Image
Image('http://www.rhsmpsychology.com/images/monocular_IV.jpg')

Out[4]:

Extending Olshausens classical SparseNet

Laurent Perrinet

2017-03-17 10:49

In a previous notebook, we tried to reproduce the learning strategy specified in the framework of the SparseNet algorithm from Bruno Olshausen. It allows to efficiently code natural image patches by constraining the code to be sparse. In particular, we saw that in order to optimize competition, it is important to control cooperation and we implemented a heuristic to just do this.
In this notebook, we provide an extension to the SparseNet algorithm. We will study how homeostasis (cooperation) may be an essential ingredient to this algorithm working on a winner-take-all basis (competition). This extension has been published as Perrinet, Neural Computation (2010) (see https://laurentperrinet.github.io/publication/perrinet-10-shl ):

@article{Perrinet10shl,
    Title = {Role of homeostasis in learning sparse representations},
    Author = {Perrinet, Laurent U.},
    Journal = {Neural Computation},
    Year = {2010},
    Doi = {10.1162/neco.2010.05-08-795},
    Keywords = {Neural population coding, Unsupervised learning, Statistics of natural images, Simple cell receptive fields, Sparse Hebbian Learning, Adaptive Matching Pursuit, Cooperative Homeostasis, Competition-Optimized Matching Pursuit},
    Month = {July},
    Number = {7},
    Url = {https://laurentperrinet.github.io/publication/perrinet-10-shl},
    Volume = {22},
}

This is joint work with Victor Boutin.

Reproducing Olshausen's classical SparseNet (part 3)

Laurent Perrinet

2017-03-16 13:51

In this notebook, we test the convergence of SparseNet as a function of different learning parameters. This shows the relative robustness of this method according to the coding parameters, but also the importance of homeostasis to obtain an efficient set of filters:

first, whatever the learning rate, the convergence is not complete without homeostasis,
second, we achieve better convergence for similar learning rates and on a certain range of learning rates for the homeostasis
third, the smoothing parameter alpha_homeo has to be properly set to achieve a good convergence.
last, this homeostatic rule works with the different variants of sparse coding.

Reproducing Olshausens classical SparseNet (part 2)

Laurent Perrinet

2017-03-15 16:01

In a previous notebook, we tried to reproduce the learning strategy specified in the framework of the SparseNet algorithm from Bruno Olshausen. It allows to efficiently code natural image patches by constraining the code to be sparse.
However, the dictionaries are qualitatively not the same as the one from the original paper, and this is certainly due to the lack of control in the competition during the learning phase.
Herein, we re-implement the cooperation mechanism in the dictionary learning routine - this will be then proposed to the main code.

This is joint work with Victor Boutin.

Reproducing Olshausen's classical SparseNet (part 1)

Laurent Perrinet

2017-03-14 13:51

This notebook tries to reproduce the learning strategy specified in the framework of the SparseNet algorithm from Bruno Olshausen. It allows to efficiently code natural image patches by constraining the code to be sparse.
the underlying machinery uses a similar dictionary learning as used in the image denoising example from sklearn and our aim here is to show that a novel ingredient is necessary to reproduce Olshausen's results.
All these code bits is regrouped in the SHL scripts repository (where you will also find some older matlab code). You may install it using

    pip install git+https://github.com/bicv/SHL_scripts

following this failed PR to sklearn that was argued in this post (and following) the goal of this notebooks is to illustrate the simpler code implemented in the SHL scripts

This is joint work with Victor Boutin.

Some basics around probabilities

Laurent Perrinet

2017-03-09 11:41

basics of probability theory¶

In the context of a course in Computational Neuroscience, I am teaching a basic introduction in Probabilities, Bayes and the Free-energy principle.

Let's learn to use probabilities in practice by generating some "synthetic data", that is by using the computer's number generator. 2017-03-13_NeuroComp_FEP

Bogacz (2017) A tutorial on free-energy

Laurent Perrinet

2017-01-15 21:41

I enjoyed reading "A tutorial on the free-energy framework for modelling perception and learning" by Rafal Bogacz, which is freely available here. In particular, the author encourages to replicate the results in the paper. He is himself giving solutions in matlab, so I had to do the same in python all within a notebook...

Laurent Perrinet

2016-11-17 09:56

Sometimes, you need to pick up the $N$-th extremal values in a mutli-dimensional matrix.

Let's suppose it is represented as a nd-array (here, I further suppose you are using the numpy library from the python language). Finding extremal values is easy with argmax, argmin or argsort but this function operated on 1d vectors... Juggling around indices is sometimes not such an easy task, but luckily, we have the unravel_index function.

For those in a hurry, one quick application is that given an np.ndarray, it's eeasy to get the index of the maximal value in that array :

In [1]:

import numpy as np
x = np.arange(2*3*4).reshape((4, 3, 2))
x = np.random.permutation(np.arange(2*3*4)).reshape((4, 3, 2))
print ('input ndarray', x)
idx = np.unravel_index(np.argmax(x.ravel()), x.shape)
print ('index of maximal value = ', idx, ' and we verify that ', x[idx], '=', x.max())

input ndarray [[[ 4  3]
  [19  7]
  [ 1 13]]

 [[ 6 12]
  [14  0]
  [16 11]]

 [[15 17]
  [20 22]
  [ 9  5]]

 [[ 8 18]
  [23 21]
  [10  2]]]
index of maximal value =  (3, 1, 0)  and we verify that  23 = 23

Let's unwrap how we found such an easy solution...

Saving and displaying movies and dynamic figures

Laurent Perrinet

2016-11-15 09:56

It is insanely useful to create movies to illustrate a talk, blog post or just to include in a notebook:

In [1]:

from IPython.display import HTML
HTML('<center><video controls autoplay loop src="../files/2016-11-15_noise.mp4" width=61.8%/></center>')

Out[1]:

For years I have used a custom made solution made around saving single frames and then calling ffmpeg to save that files to a movie file. That function (called anim_save had to be maintained accross different libraries to reflect new needs (going to WEBM and MP4 formats for instance). That made the code longer than necessary and had not its place in a scientific library.

Here, I show how to use the animation library from matplotlib to replace that

A pi-pie named Monte Carlo

Laurent Perrinet

2016-10-01 14:22

A lively community of people including students, researchers and tinkerers from Marseille (France) celebrate the so-called "π-day" on the 3rd month, 14th day of each year. A nice occasion for general talks on mathematics and society in a lively athmosphere and of course to ... a pie contest!

I participated last year (in 2016) with a pie called "Monte Carlo". Herein, I give the recipe by giving some clues on its design... This page is a notebook - meaning that you can download it and re-run the analysis I do here at home (and most importantly comment or modify it and correct potential bugs...).

Une tarte au pi nommée Monte Carlo

Laurent Perrinet

2016-10-01 14:22

Une active communauté d'étudiants, chercheurs et bidouilleurs célèbrent à Marseille la "journée π" le 3ème mois, 14ème jour de chaque année. Une occasion de rêve pour en apprendre plus sur les mathématiques et la science dans une ambiance conviviale... Mais c'est aussi l'occasion d'un concours de tartes!

J'ai eu l'opportunité d'y participer l'an dernier (soit pour l'édition 2016) avec une tarte appelée "Monte Carlo". Je vais donner ici la "recette" de ma tarte, le lien avec le nombre π et quelques digressions mathématiques (notament par rapport à la présence incongrue d'un éléphant mais aussi par rapport à la démarche scientifique)... Cette page est un "notebook" - vous pouvez donc la télécharger et relancer les analyses et figures (en utilisant python + jupyter). C'est aussi un travail non figé - prière de me suggérer des corrections!

Predictive coding of motion in an aperture

Laurent Perrinet

2016-07-16 07:58

After reading the paper Motion Direction Biases and Decoding in Human Visual Cortex by Helena X. Wang, Elisha P. Merriam, Jeremy Freeman, and David J. Heeger (The Journal of Neuroscience, 10 September 2014, 34(37): 12601-12615; doi: 10.1523/JNEUROSCI.1034-14.2014), I was interested to test the hypothesis they raise in the discussion section :

The aperture-inward bias in V1–V3 may reflect spatial interactions between visual motion signals along the path of motion (Raemaekers et al., 2009; Schellekens et al., 2013). Neural responses might have been suppressed when the stimulus could be predicted from the responses of neighboring neurons nearer the location of motion origin, a form of predictive coding (Rao and Ballard, 1999; Lee and Mumford, 2003). Under this hypothesis, spatial interactions between neurons depend on both stimulus motion direction and the neuron's relative RF locations, but the neurons themselves need not be direction selective. Perhaps consistent with this hypothesis, psychophysical sensitivity is enhanced at locations further along the path of motion than at motion origin (van Doorn and Koenderink, 1984; Verghese et al., 1999).

Concerning the origins of aperture-inward bias, I want to test an alternative possibility. In some recent modeling work:

Laurent Perrinet, Guillaume S. Masson. Motion-based prediction is sufficient to solve the aperture problem. Neural Computation, 24(10):2726--50, 2012 https://laurentperrinet.github.io/publication/perrinet-12-pred

I was surprised to observe a similar behavior: the trailing edge was exhibiting a stronger activation (i. e. higher precision revealed by a lower variance in this probabilistic model) while I would have thought intuitively the leading edge would be more informative. In retrospect, it made sense in a motion-based prediction algorithm as information from the leading edge may propagate in more directions (135° for a 45° bar) than in the trailing edge (45°, that is a factor of 3 here). While we made this prediction we did not have any evidence for it.

In this script the predictive coding is done using the MotionParticles package and for a motion texture within a disk aperture.

Static motion clouds?

Laurent Perrinet

2016-07-14 15:18

Motion Clouds were originally defined as moving textures controlled by a few parameters. The library is also capable to generate a static spatial texture. Herein I describe a solution to generate a single static frame.

compiling notebooks into a report

Laurent Perrinet

2016-06-25 07:58

For a master's project in computational neuroscience, we adopted a quite novel workflow to go all the steps from the learning of the small steps to the wrtiting of the final thesis. Though we were flexible in our method during the 6 months of this work, a simple workflow emerged that I describe here.

Compiling a set of notebook to a LaTeX document.

Compiling and using pyNN + NEST + python3

Laurent Perrinet

2016-06-01 21:21

PyNN is a neural simulation language which works well with the NEST simulator. Here I show my progress in using both with python 3 and how to show results in a notebook.

A bit more fun with gravity waves

Laurent Perrinet

2016-04-24 16:32

More fun with gravity waves¶

Motion Clouds were defined in the origin to provide a simple parameterization for textures. Thus we used a simple unimodal, normal distribution (on the log-radial frequency space to be more precise). But the larger set of Random Phase Textures may provide some interesting examples, some of them can even be fun! This is the case of this simulation of the waves you may observe on the surface on the ocean.

Main features of gravitational waves are:

longer waves travel faster (tsunami are fast and global, ripples are slow and local) - speed is linearly proportional to wavelength
phase speed (following a wave's crest) is twice as fast as group speed (following a group of waves).

More info about deep water waves : http://farside.ph.utexas.edu/teaching/336L/Fluidhtml/node122.html

IRM clouds

Laurent Perrinet

Laurent Perrinet

2016-01-19 14:03

L'installation Elasticité dynamique agit comme un filtre et génère de nouveaux espaces démultipliés, comme un empilement quasi infini d'horizons. Par principe de réflexion, la pièce absorbe l'image de l'environnement et accumule les points de vue ; le mouvement permanent requalifie continuellement ce qui est regardé et entendu.

On va maintenant utiliser des forces elastiques pour coordonner la dynamique des lames dans la trame.

bootstraping posts for élasticité

Laurent Perrinet

2016-01-18 14:03

L'installation Elasticité dynamique agit comme un filtre et génère de nouveaux espaces démultipliés, comme un empilement quasi infini d'horizons. Par principe de réflexion, la pièce absorbe l'image de l'environnement et accumule les points de vue ; le mouvement permanent requalifie continuellement ce qui est regardé et entendu.

.. media:: http://vimeo.com/150813922

Ce meta-post gère la publication sur https://laurentperrinet.github.io/sciblog.

Reproducing Olshausen's classical SparseNet (part 3)

Laurent Perrinet

2015-12-11 13:51

This is an old blog post, see the newer version in this post

2015-09-30 14:03

L'installation Elasticité dynamique agit comme un filtre et génère de nouveaux espaces démultipliés, comme un empilement quasi infini d'horizons. Par principe de réflexion, la pièce absorbe l'image de l'environnement et accumule les points de vue ; le mouvement permanent requalifie continuellement ce qui est regardé et entendu.

Ce post simule un rendu de reflection utilisant povray.

A scene with mirrors rendered with vapory (see https://laurentperrinet.github.io/sciblog/posts/2015-01-16-rendering-3d-scenes-in-python.html )

Smooth transition between MCs

Laurent Perrinet

2015-07-02 15:06

A smooth transition of MotionClouds while smoothly changing their parameters.

A simpler e-mail client: Diving into mutt

Laurent Perrinet

2015-05-29 10:55

My installation notes for mutt_+Homebrew_+gmail, based on this post by steve losh and this other post

A hitchhiker guide to Matching Pursuit

Laurent Perrinet

2015-05-22 10:59

The Matching Pursuit algorithm is popular in signal processing and applies well to digital images.

I have contributed a python implementation and we will show here how we may use that for extracting a sparse set of edges from an image.

this will be exposed in the following book chapter (see also https://laurentperrinet.github.io/publication/perrinet-15-bicv ), to be published by the end of this year:

@inbook{Perrinet15bicv,
    title = {Sparse models},
    author = {Perrinet, Laurent U.},
    booktitle = {Biologically-inspired Computer Vision},
    chapter = {13},
    editor = {Keil, Matthias and Crist\'{o}bal, Gabriel and Perrinet, Laurent U.},
    publisher = {Wiley, New-York},
    year = {2015}
}

A simple pre-processing filter for image processing

Laurent Perrinet

2015-05-21 22:01

When processing images, it is useful to avoid artifacts, in particular when you try to understand biological processes. In the past, I have used natural images (found on internet, grabbed from holiday pictures, ...) without controlling for possible problems.

In particular, digital pictures are taken on pixels which are most often placed on a rectangular grid. It means that if you rotate that image, you may lose information and distort it and thus get wrong results (even for the right algorithm!). Moreover, pictures have a border while natural scenes do not, unless you are looking at it through an aperture. Intuitively, this means that large objects would not fit on the screen and are less informative.

In computer vision, it is easier to handle these problems in Fourier space. There, an image (that we suppose square for simplicity) is transformed in a matrix of coefficients of the same size as the image. If you rotate the image, the Fourier spectrum is also rotated. But as you rotate the image, the information that was in the corners of the original spectrum may span outside the spectrum of the rotated image. Also, the information in the center of the spectrum (around low frequencies) is less relevant than the rest.

Here, we will try to keep as much information about the image as possible, while removing the artifacts related to the process of digitalizing the picture.

Extending Olshausens classical SparseNet

Laurent Perrinet

2015-05-12 10:49

This is an old blog post, see the newer version in this post

Reproducing Olshausens classical SparseNet (part 2)

Laurent Perrinet

2015-05-06 16:01

This is an old blog post, see the newer version in this post

Laurent Perrinet

2015-04-26 16:35

I was recently trying to embed a clip in a jupyter notebook using MoviePy, somthing that was working smoothly with python 2.7. After switching to python 3, this was not working anymore and left me scratcthing my head for a solution.

We live in a open-sourced world, so I filled an issue:

https://github.com/Zulko/moviepy/issues/160

elastic-force

Laurent Perrinet

2015-04-10 14:03

L'installation Elasticité dynamique agit comme un filtre et génère de nouveaux espaces démultipliés, comme un empilement quasi infini d'horizons. Par principe de réflexion, la pièce absorbe l'image de l'environnement et accumule les points de vue ; le mouvement permanent requalifie continuellement ce qui est regardé et entendu.

On va maintenant utiliser des forces elastiques pour coordonner la dynamique des lames dans la trame.

The Vancouver set

Laurent Perrinet

2015-03-26 11:06

I have heard Vancouver can get foggy and cloudy in the winter. Here, I will provide some examples of realistic simulations of it...

This stimulation was used in the following poster presented at VSS:


@article{Kreyenmeier2016,
author = {Kreyenmeier, Philipp and Fooken, Jolande and Spering, Miriam},
doi = {10.1167/16.12.457},
issn = {1534-7362},
journal = {Journal of Vision},
month = {sep},
number = {12},
pages = {457},
publisher = {The Association for Research in Vision and Ophthalmology},
title = {{Similar effects of visual context dynamics on eye and hand movements}},
url = {http://jov.arvojournals.org/article.aspx?doi=10.1167/16.12.457},
volume = {16},
year = {2016}
}

elastic-grids-of-edges

Laurent Perrinet

2015-02-26 14:03

Dans un notebook précédent, on a vu comment créer une grille hexagonale et comment l'animer.

On va maintenant utiliser MoviePy pour animer ces plots.

Saving files in HoloViews

Laurent Perrinet

2015-02-25 13:38

Here, we show how to save files from figures generated with HoloViews. This was thoroughly explained in this page.

An optic flow with Motion Clouds

Laurent Perrinet

2015-02-19 14:38

Horizon¶

Script done in collaboration with Jean Spezia.

Laurent Perrinet

2015-01-16 14:36

The above snippet shows how you can create a 3D rendered scene in a few lines of codes (from http://zulko.github.io/blog/2014/11/13/things-you-can-do-with-python-and-pov-ray/):

In [1]:

import vapory

camera = vapory.Camera( 'location', [0, 2, -3], 'look_at', [0, 1, 2] )
light = vapory.LightSource( [2, 4, -3], 'color', [1, 1, 1] )
sphere = vapory.Sphere( [0, 1, 2], 2, vapory.Texture( vapory.Pigment( 'color', [1, 0, 1] )))

scene = vapory.Scene(camera = camera , # a Camera object
                     objects = [light, sphere], # POV-Ray objects (items, lights)
                     included = ["colors.inc"]) # headers that POV-Ray may need

# passing 'ipython' as argument at the end of an IPython Notebook cell
# will display the picture in the IPython notebook.
scene.render('ipython', width=900, height=500)

Out[1]:

Here are more details...

The right imports in a notebook

Laurent Perrinet

2015-01-07 14:14

Following this post http://carreau.github.io/posts/10-No-PyLab-Thanks.ipynb.html, here is ---all in one single cell--- the bits necessary to import most useful libraries in an ipython notebook:

In [1]:

# import numpy and set the printed precision to something humans can read
import numpy as np
np.set_printoptions(precision=2, suppress=True)
# set some prefs for matplotlib
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams.update({'text.usetex': True})
fig_width_pt = 700.  # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27               # Convert pt to inches
fig_width = fig_width_pt*inches_per_pt  # width in inches
FORMATS = ['pdf', 'eps']
phi = .5*np.sqrt(5) + .5 # useful ratio for figures
# define plots to be inserted interactively
%matplotlib inline
#%config InlineBackend.figure_format='retina' # high-def PNGs, quite bad when using file versioning
%config InlineBackend.figure_format='svg'

Below, I detail some thoughts on why it is a perfect preamble for most ipython notebooks.

todo.txt

Laurent Perrinet

2014-12-29 22:05

In a previous post, I described steps to follow to live with decentralized, open-source cloud services (see here), let's focus on setting up a todo list.

Everything revolves under the http://todotxt.com/ specifications -- all in one, simple todo.txt file

Bib Cloud

Laurent Perrinet

2014-12-29 15:48

I always used a single babel.bib BibTex file to keep a record of all my readings. Great. Then I managed it using a VCS (first SVN, then git). Super great. But at the same time, resources like CiteUlike or Mendeley (but also other services like orcid) allow for cloud-like services of having this data everywhere, anytime. Super super great!

But this fails if you have no connection to internet (remote conference, ...) or more importantly if these services change their policy (mendeley to citeulike sync disappeared at a sudden). The work you provide is not yours. These are mostly commercial services while all the open-source tools are there. Most importantly, the multiplicity of tools makes it difficult to share bibliographic data easily, while if we would have a tool to translate between them this would make it possible (without changing your habits).

Can we live without Google and co?

Laurent Perrinet

2014-12-18 16:27

The tools offered by so-called "cloud services" are useful but most often rely on the capacity to lock you in. Without leaving them completely, is there an alternative?

In a previous post on how to live with open-source cloud services and to not depend on centralized private cloud services (see here), let's focus on an install on a empty android phone.

Recruiting different population ratios in V1 using orientation components: defining a protocol

Laurent Perrinet

2014-12-10 10:56

Laurent Perrinet

2014-11-22 14:24

It's easy to record tracks (for instance while running) using the "My tracks" app on android systems. What about being able to re-use them?

In :doc:2014-11-22-reading-kml-my-tracks-files-in-ipython we reviewed different approches. Let's now try to use this data.

Reading KML My Tracks files in ipython

Laurent Perrinet

2014-11-15 14:38

It's easy to record tracks (for instance while running) using the "My tracks" app on android systems. What about being able to re-use them?

Here, I am reviewing some methods to read a KML file, including fastkml, pykml, to finally opt for a custom method with the xmltodict package.

Reverse-phi and Asymmetry of ON and OFF responses

Laurent Perrinet

2014-11-10 11:03

We test Reverse-phi motion and the Asymmetry of ON and OFF responses using MotionClouds.

Recruiting different population ratios in V1 using orientation components

Laurent Perrinet

2014-11-10 10:56

A feature of MotionClouds is the ability to precisely tune the precision of information following the principal axes. One which is particularly relevant for the primary visual cortical area of primates (area V1) is to tune the orientation mean and bandwidth.

Recruiting different population ratios in V1 using orientation components: a biphoton study

Laurent Perrinet

2014-11-10 10:56

A feature of MotionClouds is the ability to precisely tune the precision of information following the principal axes. One which is particularly relevant for the primary visual cortical area of primates (area V1) is to tune the otirentation mean and bandwidth.

To install the necessary libraries, check out the documentation.

summary of the biphoton protocol¶

For the biphoton experiment:

The refresh rate of the screen is 70Hz and stimulation for 5 times 1s, which makes 350 images.
for the spatial frequency 0.125 cyc/deg is optimal (between 0.01 and 0.16).
for the temporal frequency 2 cyc/sec is optimal (between 0.8 and 4 sic/sec), we manipulate $B_V$ to get a qualitative estimate.

A bit of fun with gravity waves

Laurent Perrinet

2014-10-24 16:32

A bit of fun with gravity waves¶

Motion Clouds were defined in the origin to provide a simple parameterization for textures. Thus we used a simple unimodal, normal distribution (on the log-radial frequency space to be more precise). But the larger set of Random Phase Textures may provide some interesting examples, some of them can even be fun! This is the case of this simulation of the waves you may observe on the surface on the ocean.

In [1]:

from IPython.display import HTML
HTML('<center><video controls autoplay loop src="../files/2014-10-24_waves/waves.mp4" width=61.8%/></center>')

Out[1]:

Main features of gravitational waves are:

longer waves travel faster (tsunami are fast and global, ripples are slow and local) - speed is linearly proportional to wavelength
phase speed (following a wave's crest) is twice as fast as group speed (following a group of waves).

More info about deep water waves : http://farside.ph.utexas.edu/teaching/336L/Fluidhtml/node122.html

Statistics of time-varying images

Laurent Perrinet

2014-10-21 13:36

In this notebook, we are interested in replicating the results from Dong and Attick (1995).

Motion Plaids

Laurent Perrinet

2014-10-20 17:50

MotionPlaids : from MotionClouds components to a plaid-like stimulus¶

In [1]:

import numpy as np
import MotionClouds as mc
downscale = 2
fx, fy, ft = mc.get_grids(mc.N_X/downscale, mc.N_Y/downscale, mc.N_frame/downscale)
name = 'MotionPlaid'

mc.figpath = '../files/2014-10-20_MotionPlaids'
import os
if not(os.path.isdir(mc.figpath)): os.mkdir(mc.figpath)

Plaids are usually created by adding two moving gratings (the components) to form a new simulus (the pattern). Such stimuli are crucial to understand how for instance information about motion is represented and processed and its neural mechanisms have been extensively studied notably by Anthony Movshon at the NYU. One shortcomming is the fact that these are created by components which created interference patterns when they are added (as in a Moiré pattern). The question remains to know if everything that we know about components vs pattern processing comes from these interference patterns or reaaly by the processing of the two componenents as a whole. As such, Motion Clouds are ideal candidates because they are generated with random phases: by construction, there should not be any localized interference between components. Let's verfy that:

Defintion of parameters:

In [2]:

theta1, theta2, B_theta = np.pi/4., -np.pi/4., np.pi/32

This figure shows how one can create MotionCloud stimuli that specifically target component and pattern cell. We show in the different lines of this table respectively: Top) one motion cloud component (with a strong selectivity toward the orientation perpendicular to direction) heading in the upper diagonal Middle) a similar motion cloud component following the lower diagonal Bottom) the addition of both components: perceptually, the horizontal direction is predominant.

In [3]:

print( mc.in_show_video.__doc__)


    Columns represent isometric projections of a cube. The left column displays
    iso-surfaces of the spectral envelope by displaying enclosing volumes at 5
    different energy values with respect to the peak amplitude of the Fourier spectrum.
    The middle column shows an isometric view of the faces of the movie cube.
    The first frame of the movie lies on the x-y plane, the x-t plane lies on the
    top face and motion direction is seen as diagonal lines on this face (vertical
    motion is similarly see in the y-t face). The third column displays the actual
    movie as an animation.

    Given a name, displays the figures corresponding to the Fourier spectra, the
    stimulus cubes and movies within the notebook.

Component one:

In [4]:

diag1 = mc.envelope_gabor(fx, fy, ft, theta=theta1, V_X=np.cos(theta1), V_Y=np.sin(theta1), B_theta=B_theta)
name_ = name + '_comp1'
mc.figures(diag1, name_, seed=12234565, figpath=mc.figpath)
mc.in_show_video(name_, figpath=mc.figpath)

Component two:

In [5]:

diag2 = mc.envelope_gabor(fx, fy, ft, theta=theta2, V_X=np.cos(theta2), V_Y=np.sin(theta2), B_theta=B_theta)
name_ = name + '_comp2'
mc.figures(diag2, name_, seed=12234565, figpath=mc.figpath)
mc.in_show_video(name_, figpath=mc.figpath)

The pattern is the sum of the two components:

In [6]:

name_ = name
mc.figures(diag1 + diag2, name, seed=12234565, figpath=mc.figpath)
mc.in_show_video(name_, figpath=mc.figpath)

A grid of plaids shows switch from transparency to pattern motion¶

Script done in collaboration with Jean Spezia.

In [7]:

import os
import numpy as np
import MotionClouds as mc

name = 'MotionPlaid_9x9'
if mc.check_if_anim_exist(name):
    B_theta = np.pi/32
    N_orient = 9
    downscale = 4
    fx, fy, ft = mc.get_grids(mc.N_X//downscale, mc.N_Y//downscale, mc.N_frame//downscale)
    mov = mc.np.zeros(((N_orient*mc.N_X//downscale), (N_orient*mc.N_X//downscale), mc.N_frame//downscale))
    i = 0
    j = 0
    for theta1 in np.linspace(0, np.pi/2, N_orient)[::-1]:
        for theta2 in np.linspace(0, np.pi/2, N_orient):
            diag1 = mc.envelope_gabor(fx, fy, ft, theta=theta1, V_X=np.cos(theta1), V_Y=np.sin(theta1), B_theta=B_theta)
            diag2 = mc.envelope_gabor(fx, fy, ft, theta=theta2, V_X=np.cos(theta2), V_Y=np.sin(theta2), B_theta=B_theta)
            mov[(i)*mc.N_X//downscale:(i+1)*mc.N_X//downscale, (j)*mc.N_Y//downscale : (j+1)*mc.N_Y//downscale, :] = mc.random_cloud(diag1 + diag2, seed=1234)
            j += 1
        j = 0
        i += 1
    mc.anim_save(mc.rectif(mov, contrast=.99), os.path.join(mc.figpath, name))

In [8]:

mc.in_show_video(name, figpath=mc.figpath)

As in (Rust, 06) we show in this table the concatenation of a table of 9x9 MotionPlaids where the angle of the components vary on respectively the horizontal and vertical axes.
The diagonal from the bottom left to the top right corners show the addition of two component MotionClouds of similar direction: They are therefore also intance of the same Motion Clouds and thus consist in a single component.
As one gets further away from this diagonal, the angle between both component increases, as can be seen in the figure below. Note that first and last column are different instance of similar MotionClouds, just as first and last lines in in the table.

A parametric exploration on MotionPlaids¶

In [9]:

N_orient = 8
downscale = 2
fx, fy, ft = mc.get_grids(mc.N_X/downscale, mc.N_Y/downscale, mc.N_frame)
theta = 0
for dtheta in np.linspace(0, np.pi/2, N_orient):
    name_ = name + 'dtheta_' + str(dtheta).replace('.', '_')
    diag1 = mc.envelope_gabor(fx, fy, ft, theta=theta + dtheta, V_X=np.cos(theta + dtheta), V_Y=np.sin(theta + dtheta), B_theta=B_theta)
    diag2 = mc.envelope_gabor(fx, fy, ft, theta=theta - dtheta, V_X=np.cos(theta - dtheta), V_Y=np.sin(theta - dtheta), B_theta=B_theta)
    mc.figures(diag1 + diag2, name_, seed=12234565, figpath=mc.figpath)
    mc.in_show_video(name_, figpath=mc.figpath)

Is there any sort of utility I can use to convert an image to ASCII and then print it in my terminal? I looked for one but couldn't seem to find any.

transferring (lots of) files to a remote server + saving space

Laurent Perrinet

2014-07-07 10:31

I have recently been asked how to transfer lots of file from a backup server to a local disk. The context is that

the backup is an archive to be put on an external hard drive and then in a closet for future use,
the transfer line is not quite reliable and transfers may stop (for instance if it is mounted as a samba / windows / applefile share)

animation-in-a-notebook

Laurent Perrinet

2014-06-25 14:03

animation in a notebook¶

Dans le notebook précédent, on a vu comment créer

la grille rectangulaire puis hexagonale,
les lames autour de ces points - puis un fonction pour les représenter.

On va maintenant utiliser:

http://matplotlib.org/api/animation_api.html

http://jakevdp.github.io/blog/2012/08/18/matplotlib-animation-tutorial/

La densité est définie sur $[0, 2\pi[$ par $ f(x) = \frac {e^{\kappa cos(\theta - m)}}{2 \pi I_{0}(\kappa)} ~$
avec $~ \kappa = \frac {1}{\sigma^{2}} $

L'orientation est définie sur $[0, \pi[$ par un remplaçant de la variable $~ (\theta_d - m) ~$ par $~ 2(\theta_o - m)$

nous obtenons alors la densité pour $\theta_o$: $ f(x) = \frac{1}{2 \pi I_{0}(\kappa)} \cdot {e^{\frac{cos(2(\theta - m))}{\sigma^{2}}}}$ et $~ p = \frac {1} {2\pi \sigma_1 \sigma_2} \cdot e^{- 2 \frac{(m_2-m_1)^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}}$

car avec le changement de variable : $ \varepsilon(x) = \frac {1} {\sigma_{1} \sqrt{2\pi}} \cdot e^{\frac {-(2(x-m_{1}))^{2}} {2\sigma_{1}^{2}}} $ et $\gamma (x) = \frac {1} {\sigma_2 \sqrt{2\pi}} \cdot e^{\frac {-(2(x-m_2))^{2}} {2\sigma_{2}^{2}}} $

donc $\varepsilon(x) \cdot \gamma(x) = \frac {1}{2\pi \sigma_{1} \sigma_{2}} \cdot e^{-2(\frac {(x-m_{1})^{2}}{\sigma_{1}^{2}} + \frac{(x-m_{2})^{2}}{\sigma_{2}^{2}})} = \frac {1} {2\pi \sigma_1 \sigma_2} \cdot e^{- 2 \frac{(m_2-m_1)^{2}}{\sigma_{1}^{2} + \sigma_{2}^{2}}}$

[Stage M1] Chloé Pasturel : week 4

Laurent Perrinet

2014-06-18 11:20

simulations avec différents MCs
ring avec représentation de la vitesse angulaire

[Stage M1] Chloé Pasturel : week 3

Laurent Perrinet

2014-06-18 11:20

Travail sur les Motion Clouds et observation des différents changement de paramètre particulièrement B_sf

introduction aux motions clouds: installation et display dans un notebook
synthèse de clouds avec différents theta
synthèse avec différents B_theta (V=0)
utilisation du trick dans https://laurentperrinet.github.io/sciblog/posts/2015-07-02_smooth-transition-between-mcs.html pour créer un stimulus pour lequel l'orientation tourne de 0 à 2*pi
proposer une façon simple de passer de ce signal à une entré pour le ring (un peu de maths? un MC est une texture définie en fourier, une cellule simple fait une convolution dans l'espace, donc une multiplication dans Fourier: on pourrait avoir la sortie linéaire du neurone de V1 directement ...)

http://neuralensemble.github.io/MotionClouds/

https://neuralensemble.github.io/MotionClouds

[Stage M1] Chloé Pasturel : week 2

Laurent Perrinet

2014-06-18 11:20

Simulation d'un ring model grâce au programme brian et observation des conséquences du changement de différents paramêtres

simulation d'un ring par l'intermédiaire de brian http://briansimulator.org/
courbe de sélectivité en fonction des parametres
courbe de sélectivité à différents contrastes

[Stage M1] Chloé Pasturel : plan du stage

Laurent Perrinet

2014-06-18 11:10

Modele d'une entrée naturelle à un modèle de selectivité à l'orientation¶

Stage M1 de Chloé Pasturel, encadrée par Laurent Perrinet (INT, CNRS)

Contexte du stage¶

(importé depuis https://laurentperrinet.github.io/grant/anr-bala-v1/ )

[Stage M1] Chloé Pasturel : week 1

Laurent Perrinet

2014-06-18 11:10

définition du sujet - lecture Hansel
mise en place du python notebook
demos de brian

installation des outils¶

virtual box
VM ubuntu
<<< guest additions + clone
neurodebian

    wget -O- http://neuro.debian.net/lists/trusty.de-md.full | sudo tee /etc/apt/sources.list.d/neurodebian.sources.list
    sudo apt-key adv --recv-keys --keyserver pgp.mit.edu 2649A5A9
    sudo apt-get update

install on ubuntu

sudo apt-get install aptitude
sudo aptitude install ipython-notebook
sudo aptitude install python-matplotlib
sudo aptitude install python-brian 
sudo aptitude install texlive-full


sudo aptitude install python-pynn 
pip install --user neurotools

Batch removing files on a FTP server

Laurent Perrinet

2014-02-28 09:33

In a previous post, I have shown how to convert MoinMoin pages to this blog (Nikola angine). Let's now tidy up the place and remove obsolete pages.

Pages that were succesfully converted:

Converting moinmoin pages to nikola

Laurent Perrinet

2014-02-15 15:09

nikolFrom a previous post, we have this function to import a page from moinmoin, convert it and publish:

Converting MoinMoin pages to Nikola

Laurent Perrinet

2014-02-11 16:07

I have another wiki to take notes. These are written using the MoinMoin syntax, which is nice, but I found no converter from MoinMoin to anything compatible with Nikola.

Following this post http://carreau.github.io/posts/06-NBconvert-Doc-Draft.html , I thought I might give it a try using ipython within nikola:

In [1]:

URL = 'https://URL/cgi-bin/index.cgi/SciBlog/2005-10-29?action=print'

NOTE: this should not be tried with these URLS as they do not exist anymore...

Trying out ipython blogging

Laurent Perrinet

2014-01-14 10:37

How to easily publish an ipython notebook?¶

installing the appropriate plugins¶

Here, we will need to install additions to the nikola publishing tool: a rendering machine + a theme. Using homebrew in mcosx, this looks like:

brew install npm
npm install -g less

and

nikola install_theme zen-ipython

Once this was done, one could tune conf.py and then issue:

nikola new_post -f ipynb

Finally, one needs to build (nikola build) and publish (nikola deploy)

example¶

In [1]:

plot(rand(50))

Out[1]:

[<matplotlib.lines.Line2D at 0x10ab2d1d0>]

2013-11-12 craac

Laurent Perrinet

2012-07-09 22:28

craac 2013

Compte rendu annuel d'activité des chercheurs du CNRSAnnée 2012 - 2013

Implementing a retinotopic transform using grid_sample from pyTorch¶

SparseEdges : computing sparseness of natural images with retina-like RFs¶

Initialization¶

one example image¶

Running simulations on a set of natural images¶

Analysing results¶

name_database='serre07_distractors'¶

trying different fits¶

on the modulation function¶

on the pdf¶

TODO : do the same for experiment + '_linear'¶

DropLets¶

modeling different sparseness exponents in a natural image¶

some book keeping for the notebook¶

Poster GDR Vision¶

basics of probability theory¶

generating databases¶

More fun with gravity waves¶

Studying the role of contrast in V1 using MotionClouds¶

Horizon¶

summary of the electro-physiology protocol¶

summary of the biphoton protocol¶

A bit of fun with gravity waves¶

MotionPlaids : from MotionClouds components to a plaid-like stimulus¶

A grid of plaids shows switch from transparency to pattern motion¶

Implementing a retinotopic transform using `grid_sample` from `pyTorch`¶